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Chapter  1 
INTRODUCTION 

There  are  three  Boundary-Value  Proble;.is  (BVP)  of  poten¬ 
tial  theory  (Moritz,  1964,  p.l;  Heiskanen  and  Moritz,  1967, 
p . 36 ) : 

1.  First  BVP  (Dirichlet's  problem):  Given  an  arbitrary 
function  on  a  surface  S  ,  determine  a  function  V  which 
is  harmonic  either  inside  or  outside  S  ,  and  which 
assiunes  on  S  the  values  of  the  prescribed  function. 

2.  Second  BVP  (Neumann's  problem):  Given  the  normal  derivative 
3V/3n  of  a  function  on  a  surface  S  ,  determine  the 
function  V  which  is  harmonic  either  inside  or  outside 

S  ,  whose  normal  derivative  assumes  the  prescribed  boundary 
values  on  S  . 

3.  Third  BVP:  Determine  a  function  V  which  is  harmonic 
either  inside  or  outside  a  given  surface  S  ,  and  is 
such  that  the  linear  combination  hV  +  k3V/9n  assumes 
prescribed  values  on  S  . 

It  is  clear  that  in  all  three  problems,  the  surface 
S  is  assumed  to  be  known  (i.e.  the  locations  of  points 
on  the  surface  are  determined  in  a  specified  coordinate 
system),  and  only  the  harmonic  function  V  is  to  be  determined 
from  the  various  data  on  S  .  The  geodetic  BVP  may  be  defined 
as  the  determination  of  the  physical  surface  of  the  earth, 
if  the  gravity  vector  g  and  the  gravity  potential  W  are 
given  on  it  (Moritz,  1965,  p.36).  This  problem  calls  for 
the  determination  of  the  surface  of: 

(a)  the  geoid  (Stokes'  problem)  using  gravity  data  reduced 
from  the  physical  surface  of  the  earth  to  the  surface 
of  the  geoid,  or 

(b)  the  telluroid  (Molodensky ' s  problem)  from  gravity  data 
given  on  the  physical  surface  of  the  earth. 

Since  the  surface  (geoid  or  telluroid)  is  the  unknown, 
the  geodetic  BVP  is  not  directly  related  to  any  of  the  three 
problems  of  potential  theory  mentioned  above.  What  makes 
the  solution  of  the  geodetic  BVP  possible  as  a  3rd  BVP, 
is  the  fact  that  certain  restrictions  are  enforced  (Moritz, 

1964,  p.2),  namely: 

(a)  the  geoid  is  an  equipotential  surface  of  the  gravity 
field  of  the  earth  (Stokes'  approach),  and 

(b)  in  addition  to  the  gravity  g  ,  the  potential  difference 
AW  ■  W  -  Wo  at  every  surface  point  is  determined  from 

1. 
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levelling  combined  with  gravity  observations  (Molodensky 's 
approach),  Wo  being  the  potential  of  the  geoid. 

The  definitions  above  will  help  us  later  in  determining 
the  kind  of  BVP  we  need  to  solve  in  the  present  study .  Let 
us  now  state  the  problem. 

Certain  geodetic  applications,  such  as  aircraft  inertial 
navigation  systems,  and  the  computation  of  rocket  trajectories, 
or  components  of  the  gravity  disturbance  vector,  require 
the  calculation  of  the  gravity  vector  components  at  a  certain 
altitude  in  space,  from  data  collected  on  the  physical  surface 
of  the  earth.  Zero-order  solutions  to  this  problem  (i.e. 
solutions  ignoring  the  topography  of  the  earth  and/or 
the  surface  inclinations)  can  be  found  in  the  literature 
(cf.  Hirvonen  arjd  Moritz,  1963;  Heiskanen  and  Moritz,  1967; 
Mueller,  1966).  These  methods  are: 

1.  The  Direct  Integration  method,  which  is  based  on  the 
generalized  Stokes’  formula; 

2.  The  Coating  method  (see  also  in  (Orlin,  1959)); 

3.  The  Upward  Continuation  method,  which  is  based  on  the 
Poisson's  integral  equation  (Shebalin,  1979). 

Apart  from  these  three  methods,  which  are  considered 
as  t'.ie  classical  approaches,  there  are  also  other  techniques. 
However  not  all  of  them  yield  the  three  components  of  the 
gravity  disturbance  vector.  Some  of  them  are  applied  for 
the  computation  of  the  gravity  anomaly  (Ag)  only,  and  some 
others  for  the  computaion  of  the  disturbing  potential  (T) 
alone.  These  methods  are: 

4.  The  Stokes'  formula  in  combination  with  the  use  of  Poisson's 
integral  theorem  (L.  de  Witte,  1969); 

5.  The  Green's  third  identity  (Moritz,  1965;  B.  Witte,  1969; 
Koch,  196 7-b,  1968-a); 

6.  A  Bjerhammar-type  of  solution,  based  on  the  reduction 
of  gravity  data  on  a  defined  sphere  internal  to  the 
surface  of  the  earth,  and  then  application  of  the  Pois¬ 
son's  integral  for  the  upward  continuation  of  the  gravity 
data  to  a  space  point  (Bjerhammar,  1978;  Moritz,  1965, 
p.54;  Sjbberg,  1978); 

7.  Polynomical  modeling  of  the  anomalous  gravity  field, 
for  the  upward  continuation  of  gravity  data  (Paul  and 
Nagy,  1972); 

8.  Series  expansion  of  the  external  anomalous  potential 
(Petrovskaya,  1979); 

9.  The  Finite  Element  method  (Richardson  and  Hopkins,  1978; 
Junkins  and  Saunders,  1977); 
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10.  Geopotential  Modeling  by  the  so-called  point-masses 
technique  (Needham,  1970): 

11.  The  method  of  Least-Squares  Collocation,  or  a  combina¬ 
tion  of  Collocation  with  integral  equations  (Lachapelle, 
1977;  Tscherning  and  Forsberg,  1978;  Forsberg  and 
Tscherning,  1980): 

12.  The  Initial-Value  method  (Nakiboglou  and  Lim  1979). 

For  the  computation  of  the  gravity  ve' tor  g  at  a 
space  point  P  ,  we  first  need  to  evaluate  the  components 
of  the  gravity  disturbance  vector  t  at  this  point  (Heiskanen 
and  Moritz  1967,  p.227).  Since  o  =  grad  T  (more  details 
will  follow  in  the  next  chapter),  the  problem  is  equivalent 
to  the  estimation  of  the  partial  derivatives  of  the  disturbing 
potential  T  at  the  same  space  point  P  .  A  first-order 
solution,  linear  in  elevation  (h)  of  the  earth's  topography 
is  discussed  in  (Moritz,  1965.  part  2),  but  only  the  expression 
for  the  gravity  anomaly  at  P  is  derived  following  a  planar 
approximat ion . 

In  this  study  we  develop  a  method  for  the  computation 
of  the  three  components  of  the  gravity  disturbance  vector 
at  space  points,  using  data  on  the  physical  surface  of  the 
earth.  The  topography  of  the  earth  is  not  neglected,  and 
the  need  to  know  the  deflection  of  the  vertical  is  avoided. 

Only  surface  gravimetry,  elevation  data,  and  (possibly) 
a  geopotential  model  will  be  used.  Green's  third  identity 
is  the  starting  point  for  this  approach  (Chapters  2,  3, 
and  4). 

The  development  of  another  method,  based  on  the  Bjerham- 
mar's  discrete  (Dirac)  approach  for  the  upward  continuation 
of  gravity  data  is  also  investigated  (Chapters  5,  and  6). 

The  two  methods  are  then  compared  to  each  other. 

Finally,  the  applicability  of  least-squares  collocation 
is  examined,  and  the  results  of  the  three  approaches  are 
compared  with  respect  to  factors  such  as  the  inclination 
of  the  terrain,  the  altitude  of  the  space  point,  the  density 
of  the  surface  data,  etc.  (Chapter  7). 

In  order  to  avoid  the  random  and  the  systematic  errors 
which  exist  in  real  data,  we  decided  to  test  the  accuracy 
of  the  three  approaches  through  a  simulation  study,  by  substi¬ 
tuting  for  the  topography  of  the  earth  a  number  of  simple 
terrain  models  (Chapter  3).  The  exact  vectors  computed 
from  these  models  can  then  be  compared  to  the  computed  ones 
from  the  three  approaches,  using  the  synthetic  data  on  the 
model's  surface.  Such  simulation  with  terrain  models,  but 
for  different  applications  have  been  made  in  the  past  too. 

For  example  the  effect  of  the  topography  on  the  external 
gravity  anomalies  is  discussed  by  Moritz  (1965),  and  the 
computation  of  the  first  derivatives  of  the  disturbing 


potential  on  the  earth’s  surface  by  Green's  formula  is  des¬ 
cribed  by  Koch  (1967-b,  1968-a). 


Due  to  the  fact  that  throughout  this  study  the  surface 
of  the  earth  is  assumed  to  be  known,  the  problem  of  the 
determination  of  the  gravity  disturbance  bector  is  closer 
to  the  BVP  as  defined  in  potential  theory,  than  to  the  geo¬ 
detic  BVP.  Koch  and  Pope  (1972)  called  this  problem  "the 
geodetic  BVP  using  the  known  surface  of  the  earth",  and 
they  have  stated  theorems  on  its  uniqueness  and  existence. 
The  three  methods  which  will  be  described  in  the  following 
chapters,  require  as  data  (in  addition  to  the  information 
which  defines  the  surface),  surface  gravity  anomalies  (^g), 
or  disturbing  potentials  (T)  and  gravity  disturbances  (6g) 
as  boundary  values.  These  three  quantities  are  related 
to  each  other  through  the  fundamental  equation  of  physical 
geodesy  (Heiskanen  and  Moritz,  1967,  p.88):  ‘ 


Ag - 


6g  - 


If  the  boundary  values  are  gravity  anomalies  (linear 
combination  of  T  and  ^T/^r),  our  problem  is  equivalent 
to  the  3rd  BVP  above,  since  the  components  of  *  grad  T 
are  the  partial  derivatives  of  T  .  However,  if  the  bound¬ 
ary  values  are  gravity  disturbances  and  disturbing  poten¬ 
tials  (not  their  linear  combination),  our  problem  is  none 
of  the  three  BVP  of  potential  theory  as  they  are  defined 
in  the  beginning  of  this  chapter. 

The  Direct  Integration,  the  Coating,  and  the  Upward 
Continuation  methods  which  were  mentioned  above,  were  refer¬ 
red  to  as  the  classical  approaches,  because  the  topography 
of  the  earth  is  neglected,  and  the  data  required  for  their 
application  must  be  reduced  on  the  surface  of  the  sphere 
which  approximates  the  geoid.  These  classical  approaches 
are  extensively  discussed  in  (Heiskanen  and  Moritz,  1967, 
chapter  6).  The  Green’s,  the  Dirac,  and  the  collocation 
approaches  which  will  be  discussed  in  this  paper,  will 
be  referred  to  as  the  Improved  techniques,  since  the  earth's 
topography  is  not  neglected,  and  the  required  data  is  boundary 
values  on  the  earth’s  surface,  with  no  need  to  apply  gravity 
reductions  from  this  surface  to  the  geoid. 

In  chapter  8,  these  improved  techniques  are  compared 
to  the  classical  approach  (the  Direct  Integration  method). 

It  will  be  shown  that  the  results  from  the  improved  techniques 
agree  to  the  exact  values  (computed  from  the  model)  better 
than  the  results  from  the  classical  approach  do.  This 
clearly  indicates  that  these  techniques  offer  an  improved 
solution  for  the  computation  of  the  gravity  vector  in  space, 
without  neglecting  the  topography  of  the  earth. 


Chapter  2. 

THE  APPLICATION  OF  GREEN'S  THIRD  IDENTITY 


2.1  Introduction 


If  W  is  the  gravity  potential  of  the  earth,  and  U  is 
is  the  normal  gravity  potential  of  an  equipotential  ellip¬ 
soid,  then  the  earth's  gravity  vector,  and  the  normal  gravity 
vector  are  defined  as  (Heiskanen  and  Moritz,  1967,  p.85): 

t  3w  aw  - 

g  -  grad  W  -  ) 

Y  -  grad  U  =  (-^  ,  )  (2.1) 

(X,  Y,  Z)  is  a  cartesian  geocentric  coordinate  system, 
whose  origin  is  at  the  center  of  mass  of  the  earth;  the 
Z-axis  coincides  with  the  earth's  mean  rotational  axis, 
and  the  X,  Y,  axes  form  a  right-handed  system  with  the 
Z-axis,  such  that  the  X-axis  lies  in  the  Greenwich  meridian 
plane . 


The  gravity  disturbance  vector  at  any  point  P  is 
defined  as  the  difference  between  g  and  y  - 


Let  m  be  the  unit  vector  normal  to  the  equipotential 
surface  W  =  Wp  of  the  earth's  gravity  field  at  P  ,  and  I 
be  the  unit  vector  normal  to  the  equipotential  surface 
U  =  Up  of  the  normal  gravity  field  at  P  .  Then 

^  =  g  -  Y  =  grad  W  -  grad  U  =  grad  (W  -  U) 

9t  9t  9t 

=  grad  T  -  )  (2.3) 

T  being  the  disturbing  potential  at  P  .  I^  we  consider 
only  the  magnitudes  of  the  vectors  g  and  y  »  then  the 
gravity  disturbance  at  P  is  : 

5gp  -  gp  -  Yp  =  -<-§nr"  TT^P"  “^^P"  -^“gF^P 

(2.4) 


This  equation  shows  that  if  we  ignore  the  difference 
in  the  direction  of  tfi  and  ^  ,  or  in  other  words  the  tilt 
of  the  surfaces  Wp  and  Up  with  respect  to  each  other, 
then  the  normal  component  of  the  gravity  disturbance  vector 
is  identical  to  the  difference  gp  “  Yp  •  The  error  of 
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such  an  approximation  is  of  the  order  of  the  square  of  the 
deflection  of  the  vertical  (Molodensky,  et  al.,  1962,  p.79). 


The  gravity  vector  g  can  be  computed  at  any  space 
point  P  ,  provided  that  we  are  able  to  estimate  the  gravity 
disturbance  vector  t  ,  because  the  computation  of  the 
normal  gravity  vector  t  does  not  present  any  practical 
or  theoretical  difficulty  (Heiskanen  and  Moritz, 1967,  chapter 
6),  provided  that  the  altitude  of  P  is  known.  Therefore, 
if  we  evaluate  the  three  components: 

ST  3T  ST 

’  “7T  *  TT 

of  the  gravity  disturbance  vector  ^  ,  the  gravity  vector  g 
will  then  be: 


3T  3T  ST  .  .  .  SU  SU  3U  . 

''TT  ’  TT*  »  TT*  *  TT  ’  (2.5) 


2.2  Green's  Third  Identity 

If  u  denotes  the  exterior  space  of  an  arbitrary 
surface  S  ,  then.  Green’s  third  identity  is  an  integral 
equation  of  the  form  (Heiskanen  and  Moritz,  1967,  p.l2): 

///  ^  AV  du  *  -pV  -  //  -  V  ^  (i)}  dS  (2.6) 

u  S 


where ; 


V  a  continuous  and  finite  function  in  the  space 
outside  S  which  vanishes  at  infinity. 

u  the  exterior  space  of  the  surface  S  . 

n  outer  normal  to  the  surface  S  . 


the  distance  from  the  space  point  to  the  surface 
element  dS  . 


i4iT  , 
3=<2TT  , 

(o  , 


4^  ,  if  the  identity  is  applied  to  a  point  outside  S. 

2ifT  . .  . on  S . 

0  . inside  S. 


If  V  happens  to  be  a  harmonic  function  in  the  space 
outside  S  ,  by  definition  it  satisfies  Laplace's  equation 
(AV  »  0),  and  therefore,  the  left-hand  side  of  (2.6)  is 
identical  to  zero.  Furthermore,  if  the  mass  of  the  atmosphere 
is  neglected,  the  disturbing  potential  T  »  W-U  is  a 
harmonic  function  (Heiskanen  and  Moritz,  1967,  p.86): 


AT 


S^T  ^  S^T  ^ 

“TT  -TP  * 


(2.7) 


The  partial  derivatives  of  T  in  a  cartesian  coordinate 
system  are  harmonic  functions  too,  because  they  satisfy 


Laplace's  equation: 
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3^  .  3T  .  .  3237..  32  3t  . 


(2.8) 


and  similarly  for  and  — ^  .  If  we  now  apply  Green's 

identity  to  the  disturbing  potential  T  for  the  exterior 
of  the  earth's  surface  (p  =  4tt),  the  following  very  impor¬ 
tant  equation  of  physical  geodesy  is  obtained  (ibid.,  p.l2): 


T 


P 


//  fT 

S 


(2.9) 


S  being  the  physical  surface  of  the  earth,  T  the  disturb¬ 
ing  potential  at  the  surface  element  dS  ,  is  the  dis¬ 

tance  from  the  space  point  P  to  dS  ,  and  Tp  is  the 
disturbing  potential  at  P  .  Green's  third  identity  can 
also  be  applied  to  the  derivatives: 


3  T  3T  3T 
*  “7T  ' 


of  the  disturbing  potential  in  a  cartesian  coordinate  system 


(-^)p 


,  3T  .  1 

,  3T  .  1 

-Tt 


jj  I  3T  3  ,  1,  1 


"IJT  "IF  ■  T  “If  ^ 


U  ‘“I?  IF  '  T  -TH-  ^TT 


)  }  dS 


(2.10) 


f  f  J  3T  3  X  1 \  1  3  f  3T  .1  . _ 


An  initial  attempt  to  use  these  three  equations  for 
the  computation  of  the  gravity  disturbance  vector,  resulted 
in  very  complicated  expressions  the  practical  use  of  which 
was  questionable.  More  specifically,  the  deflections  of 
the  vertical  on  the  surface  S  were  needed,  as  well  as 
their  partial  derivatives  in  a  local  cartesian  coordinate 
system.  Therefore,  we  decided  to  proceed  by  differentiating 
(2.9)  directly  in  a  cartesian  coordinate  system. 

Instead  of  using  (2.9)  in  its  present  form,  certain 
authors  have  used  another  form  (Molodensky,  et  al.,  1962, 
p.45;  Koch,  1967-b,  p.29;  Koch,  1968-a,  p.ll): 

■^P  “  4f  -  T  -|f^  (2.9') 
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where  T(Po )  is  the  disturbing  potential  at  the  projection 
Pa  of  P  onto  the  surface  S  .  As  Uolodensky  (et  al., 
ibid)  explained,  this  second  form  of  (2.9)  is  an  artifact 
that  does  not  affect  the  rigour  of  the  equation.  T(Po ) 
is  a  constant,  and  the  integral 

hr  ds  »  0 

is  equal  to  zero,  if  we  integrate  over  the  whole  surface. 
However,  if  the  integration  is  taken  over  a  limited  area 
around  Po  -  say  within  a  5®  or  10®  cap  -  ,  the  term 

-  !I  T  <'“•>  lTr<f>  “s 

is  not  equal  to  zero,  and  therefore  will  cause  errors  in 
the  computation  of  Tp  .  In  our  derivation  for  the  com¬ 
ponents  of  the  gravity  disturbance  vector  in  space,  we 
decided  to  use  the  original  form  (equation  2.9),  because 
as  we  will  se  later,  a  limited  amount  of  data  is  used, 
and  not  a  global  data  set. 

Let  us  now  derive  the  equations  for  the  derivatives 
of  T  under  a  spherical  approximation  (Moritz,  1966,  p.25). 
This  does  not  mean  that  the  topography  of  the  earth  is 
not  considered,  or  neglected.  It  means  that  the  given 
heights  above  a  reference  ellipsoid,  are  taken  as  heights 
above  a  mean-earth  sphere  of  radius  R  ,  by  neglecting 
the  flattening  of  the  ellipsoid.  The  effect  of  this  spher¬ 
ical  approximation  on  the  disturbing  potential  is  estimated 
to  be  of  the  order  of  0.003  T  (Moritz,  1980,  p.l5).  Let 
Q  be  an  arbitrary  point  on  the  earth's  surface  S  ,  and 
P  the  point  in  space  where  the  derivatives  of  T  are  to 
be  computed,  (see  figure  2.1).  Let  us  also  define  the 
following  two  coordinate  systems: 

a.  A  right-handed  cartesian  coordinate  system  (X,  ?,  2) 
located  at  the  space  point  P  .  The  Z-axis  is  in  the  direc¬ 
tion  from  the  geocenter  O  to  P  ,  the  X-axis  points  to 
the  north,  and  the  ?-axis  points  to  the  west. 

b.  A  right-handed  cartesian  coordinate  system  (x,y,z) 
located  at  the  variable  surface  point  Q  (topocentric 
system).  The  z-axis  points  up,  the  x-axis  to  the  north, 
and  the  y-axis  to  the  west. 

If  we  differentiate  (2.9)  with  respect  to  X,  ?,  Z, 
we  obtain: 
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TT  Jj 

r  .  1  if 

^“5T^P'  TT  y 

^T7^p'  tt  y 


"It 

"4  ■  "Iy^T^ 

^■^“12  ■  "I^^T^ 


3T 

~5T 

3T 

"TT 

3T 

"TT 


}  dS 

Ids  (2.11) 

}  dS 


As  it  will  be  shown  in  the  next  section,  the  term 
3T/3n  is  a  function  of  the  deflection  of  the  vertical 
components  (Cil).  Hence,  (2.11)  have  the  same  difficulties 
as  (2.10)  if  applied.  For  this  reason,  equation  (2.9) 
for  T  will  be  modified  following  the  Uolodensky's  tech¬ 
nique  (Molodensky,  et  al.,  1962),  which  has  also  been  used 
in  (Moritz,  1965,  1966).  The  theory  behind  these  develop¬ 
ments  is  outlined  in  the  next  two  sections. 


2.3  Moritz's  approach 

For  a  differentiable  scalar  function  F  ,  the  gradient 
of  F  is  written  in  terms  of  its  components  as; 


grad  F  »  VF 


,  3F  3F  3F 

^"TT  '  —  '  "TT ^ 


(2.12) 


where  (x,y,z)  is  a  coordinate  system  with  arbitrary  loca¬ 
tion  and  orientation.  The  component  of  grad  F  in  the  direc¬ 
tion  of  a  unit  vector  H  is  the  dot  product: 


“It  “  "  •  ^ 


(2.13) 


Physically,  this  product  is  the  rate  of  change  of  the  scalar 
F  at  the  point  (x,y,z)  in  the  direction  of  the  vector 
n  =  (ni,  n2,  03).  Directly  from  (2.13)  we  get: 


3F 

3n 


ni 


n2 


ns 


3F 

3z 


(2.14) 


If  F  is  the  disturbing  potential  T  ,  the  term 
3F/3n  is  actually  the  term  appearing  in  (2.11).  In  this 
section,  we  follow  Moritz's  approach  for  its  evaluation 
(Moritz,  1964).  Let  (x,y,z)  be  the  local  cartesian  coor¬ 
dinate  system,  as  defined  in  the  previous  section,  located 
at  the  arbitrary  surface  point  Q  .  In  this  system,  the 
earth's  topography  can  be  represented  as 


S(x,  y,  z)  -  0 


(2.15) 


or  by  the  explicit  form  z  =  h  (x,y),  which  yields: 
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S(x  ,  y  ,  h  (x,  y) )  =  0 


(2.16) 


Differentiating  (2.16)  with  respect  to  x  and  y  we  obtain: 


dS 

3s 

3s  3h 

dx 

3x 

flS 

3S 

^  3S  3h 

-ay 

3y 

(2.17) 


If  II  is  the  unit  vector  perpendicular  to  the  surface 
S(x,y,z)  =0  at  a  point  (x,y,z),  its  components  are  propor¬ 
tional  to  the  components  of  the  vector  grad  S  : 


grad  S  =  ( 


3  S 


3x 

3S 


3  z 
3S 


3S 

3y 

3h 
3x  ’ 


3S 

3z 


) 


.  3h 

Therefore,  the  unit  vector 


3S 

3z 

3h 

-W 


3h 

3y  ’ 


3S 

3z 


) 


,  1) 

can  be  written  as; 


(2.18) 


n  = 


(--^  , 


3h 

'W'  ’ 
- 5? 


1) 


(2.19) 


If  8  is  the  angle  of  maximum  inclination  at  a  parti¬ 
cular  surface  point  Q  ,  then 

cos  6  “  a,  •  n  ^2.20 

where  n  is  the  unit  vector  normal  to  S  at  Q  ,  and 
ifo  is  the  unit  vector  along  the  z-axis; 


=  (0  ,  0  ,  1) 

n  *  (ni  ,  n2  ,  n3)  »  (ni  ,  n2 


(2.21) 


1 


i2.)2+ 
3y  ^ 


(2.22) 


Let 
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-iil 

ax 


tan 


and 


ay 


tan  Ty 


(2.23) 


Here  "fx  and  '^y  are  the  components  of  the  topographic 
inclination  in  the  x  and  y  directions  respectively. 
We  obtain  from  (2.20): 


cos  e  =  n 3  *  cos  g  -  — —  -■ 

’^tan^  Tjj  +  tan^  Xy  +  1 

and  from  well-known  trigonometric  identities: 


(2.24) 


tan^S 


tan^  +  tan^  Xy 


(2.25) 


Therefore,  the  unit  vector  n  in  (2.19)  can  be  expres¬ 
sed  as: 


n  =  cos  e  (-  tan  ,  -tan  Xy  ,  1) 


and  finally,  (2.14)  yields; 

3F  _  raF  /  3F  ^  t 


3  F 

“ajT  tan  Xy)}  cosS 


(2.26) 


(2.27) 


Equation  (2.27)  was  first  derived  by  Moritz  (1964, 
p.20),  and  it  was  applied  to  the  disturbing  potential 
T(®F)  ,  in  which  case: 


aT 

*  -yC 

(2.28-a) 

aT 

=  1-  Yn 

(2.28-b) 

aT 

(2.28-c) 

The  plus  sign  in  (2,28-b)  is  due  to  the  fact  that  in  the 
present  study  the  local  (x,y,z)  -  system  is  a  right-handed 
cartesian  system,  with  the  y-axis  pointing  to  the  west , 
instead  to  the  east  as  it  is  usually  done  in  the  literature 
(cf.  Heiskanen  and  Moritz,  1967,  p.ll2;  Moritz,  1980, 
p.l4).  5  ,  and  n  are  the  components  of  the  deflection 

of  the  vertical,  and  5g  is  the  normal  component  of  the 
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gravity  disturbance  vector.  The  normal  gravity  at  Q  is 
denoted  by  y  •  Equations  (2.28)  and  (2.27)  yield: 


-[5g  +  y(-  Ctan  n  tan  "'’y)]  cos  B 


(2.29) 


As  we  mentioned  below  equation  (2.11),  oy  substituting 

(2.29)  into  (2.11),  the  unknown  deflection  jf  the  vertical 
components  appear  in  the  integrals.  A  technique  which 
is  independent  of  these  components  has  been  developed  by 
Molodensky,  et  al.  (1962),  and  it  is  outlined  in  the  next 
section . 


2.4  Molodensky 's  Approach 

The  basic  idea  behind  this  development  is  to  modify 
the  integral  equation  (2.9)  for  the  disturbing  potential 
in  such  a  way,  that  the  term  ()T/9n  does  not  appear  any 
more  as  a  function  of  the  components  of  the  deflection 
of  the  verti  al. 


Let  us  denote  by  3T/8x  ,  3T/9y  ,  the  derivatives 

of  T  along  the  horizontal  plane  (x,y)  at  the  surface 
point  Q  ,  and  by  32T/9x  ,  ^zT/ay  ,  the  derivatives  of 
T  along  the  surface  S  (cf.  Molodensky,  et  al.,  1962, 
p.84;  Moritz,  1964,  p.21).  These  two  sets  of  partial 
derivatives  are  related  to  each  other  through  the  following 
equations: 


iil 

3x 

3T 

3x 

JJL  + 

3X 

3T 

3y 

3y  . 

3T 

dZ 

3z 

3x 

iil- 

3T 

3x  ^ 

3T 

3y  + 

3T 

3z 

3y 

3x 

3y 

3y 

3y 

3z 

3y 

or 


i,  if 

combined  with  (2. 

23) 

and  (2.: 

SzT 

3T  . 

3T 

3z 

3T 

3x 

3x 

3z 

3x 

3x 

3  ,T 

3T  . 

3T 

3z 

3T 

3y 

3z 

3y 

3y 

3T 

3x 

3x 

6g  tan 

^x 

3T 

3y 

»  hL  + 

3y 

6g  tan 

-  6g  tan  T, 


-  6g  tan  T 


(2.30) 


(2.31) 


We  can  not  substitute  (2.31)  and  (2.24)  into  (2.27) 
for  T  =  F  ,  in  order  to  obtain; 

-if  -  -  "x  *  ^  (2.32) 

This  leads  to  the  following  equation  (see  also  in  Moritz, 
1964,  p.22); 


-H  -  -  S^B-  ® 


(2.33) 


where 

D(T.h)  =  ^tan  V 

is  a  special  case  of  the  operator: 

D(U  V)  =  -^  +  ii£  ill 


(2.34) 


(2.35) 


Instead  of  substituting  (2.33)  into  (2.11)  as  we  did 
with  the  Moritz's  approach,  we  will  proceed  as  follows. 
First,  (2.33)  is  multiplied  by  1/1  to  obtain; 


—  =  -(-iS.  +  C2s_B  g  fT  h)) 

5  an  ^ilcosB  ^  (i.b)) 


(2.36) 


and  from  the  definition  of  the  operator  D  (Moritz,  1966, 
p.l8)  we  get: 


D  (T,h)  =  C  (-|-,h)  -  T  B  (4-,h) 


(2.37) 


If  now  (2.36)  is  substituted  into  the  integral  equation 
for  the  disturbing  potential,  (2.9)  yields; 

"P  =  i  IBoif-  ®  ^  (2.38) 

and  then  (2.37)  into  (2.38)  results  in: 

'^P  “  £^se'^  [D(-^,h)-TT5(-^  ,h)]}dS 

(2.39) 

Our  goal  is  to  differentiate  (2.39)  with  respect  to 
X,  ?,  and  Z  ,  but  it  is  still  necessary  to  simplify  the 
expression  for  the  operator  B  in  terms  of  known  quantities. 
We  start  from  Molodensky's  identity  (Molodensky,et  al.,  1962, 
p.85) 
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//  D(U,V)  cosB  dS  *  -/J  UA2V  cosB  dS 
S  S  (2.40) 


which  for  U  =  T/^  and  V  =  h  yields: 


//d  ,h)  cosB  dS  =  -  //^Aah  cosB  dS 

g  x.  g  *- 

where  A2  is  another  operator  defined  as  (Molodensky,  et 
al. ,  1962,  p.85) : 


A2h 


1-32  /hoh2^hv  .  32  /hohi32hv, 

hrh-2  f  “TTi  ^“ET  ^ 


(2.42) 


Note  that  in  the  last  equation  hi  ,  h2  1  and  ho  are 
functions  of  the  coordinates  Qi  »  q2  »  Qo  of  an 

orthogonal  curvilinear  system.  If  q,  and  q^  are  identi¬ 
fied  as  the  geocentric  latitude  and  longitude  respectively, 
and  q  (,  as  the  elevation,  then  a  line  element  in  this 
system  can  be  expressed  as: 

ds^  =  h§  dh^  +  hi  d^^  +  hj  dX^ 


where  (see  in  Molodensky,  et  al.,  ibid,  p.86;  Moritz, 
1966,  p.24): 


ho  =  1 
hj  =  r 
h  2  =  r  cos  (J) 
and  r  =  R  +  h 


With  this  notation,  (2.42)  now  becomes: 


A2h 


1 

r'fcos^ 


[-sin<|» 


7  SzE 


3^ 


+  cos$ 


3 |h  ^  1 

3$i  cos?! 


^h  , 


(2.43) 


If  the  first  and  second-order  derivatives  of  the  eleva¬ 
tion  at  any  surface  point  Q  are  known,  it  is  possible 
to  use  (2.39)  for  the  disturbing  potential,  or  to  differen¬ 
tiate  it  for  the  components  of  the  gravity  disturbance, 
as  it  will  be  done  in  the  next  section.  As  we  will  see, 
the  expressions  for  these  components,  not  only  require 
knowledge  of  the  elevations  h  over  the  whole  surface, 
but  they  also  require  the  slopes  tan  and  tan  ty  to 
be  known.  In  addition,  the  operator  (2.43)  requires  the 
second-order  derivatives  of  the  elevation  to  be  known  over 
the  whole  surface  S  .  At  this  point  we  will  make  the 
assumption  that  the  earth's  surface  is  a  smooth  surface. 
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whose  second-order  derivatives  (of  elevation)  are  equal 
to  zero.  This  assumption  of  a  smooth  surface  is  also  a 
condition  which  is  required  to  fulfill  Lyapunov's  conditions 
(Molodensky  et  al.,  1962,  p.83).  Under  this  assumption, 
(2.43)  becomes: 

A2h  =  -  ^  tan$  ^  ^2.44) 


For  any  function  F  defined  on  the  surface  S  ,  such 
as  the  elevation  h  ,  or  the  reciprocal  of  the  distance 
1/2.  ,  we  have  (Molodensky  et  al.,  1962,  p.84;  Moritz, 
1966,  p.22): 


32F  ^  ^ 

3qi  5qi 

3  2F  ^  3F 

'"Tq'a  3q2 


Therefore,  (2.44)  becomes: 

.  _  1  3  h 

A2h  *  -  ^  tan$ 


(2.45) 


In  (2.43)  there  is  no  first-order  derivative  of  the 
elevation  with  respect  to  longitude  X  ,  and  this  is  why 
this  term  is  missing  from  (2.45)  too.  If  we  evaluate  the 
partial  derivative  of  h  with  respect  to  the  latitude  $  , 
we  obtain: 

3h  _  3x  3h  3y 

“5T  Tx  Ti"  “Jy*  T? 


The  transformation  from  one  cartesian  coordinate  system 
to  another,  is  usually  done  by  applying  a  series  of  roata- 
tions  and  an  origin  shift  to  the  coordinate  system  being 
transformed.  In  the  present  report  we  will  use  the  following 
three  matrices  to  denote  a  positive  rotation  0  around 
each  one  of  the  three  axes: 


Ri  (0) 


10  0 

0  cos  6  sin  6 

0  -sin  6  cos  9 


] 

cos9 

0 

-sin© 

R2 

(9)  = 

0 

1 

0 

1 

sin© 

0 

cos© 

1 

cos© 

sin© 

0 

R3 

(9)  = 

-sin© 

sin© 

0 

0 

0 

1 

starting  from  the  transformation  equation  between  the 
topocentric  (x,y,z)  and  rtie  geocentric  (X,Y,Z)  systems: 


X 

X 

y 

=  Ri  [-(90»-  5q)]  Rj  [-(180°-  Xq)] 

Y 

Z 

it  is  very  easy  to  show  that: 

-  '•q  -  R  * 

» 

and  hence 


9h  _  9h  9x 
Btji  5x 


r  tan  t 


Therefore,  (2.45)  can  now  be  written  as: 

Azh  =  -  tan  t 

r  X 


(2.46) 


After  the  above  first-order  approximation  has  been 
obtained  for  Azh  ,  it  remains  to  derive  the  expression 
for  the  disturbing  potential  Tp  ,  Substituting  (2.41) 
into  (2.39)  we  obtain: 

•Tp  =  ^  -5n<T>  <=0=6 


Next,  we  apply  equation  (2.33)  with  1/1  instead  of  T 
(the  proof  is  completely  analogous)  to  obtain: 
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*  c^B  ^  ® 

and  then,  (2.47)  with  (2.48)  yield: 

"^P  “  Tit  ^4*  cosB 

■'■  5^  ^  5(-ph)- Azh  }  ds 

^  (2.49) 


As  we  have  already  mentioned,  1/S,  is  a  surface  function 
and  hence  (2.35)  yields: 


D 


isi 

3x 


tan 


+ 


ach 


tan  T 


y 


(2.50) 


As  a  final  step,  substitute  (2.50)  and  (2.46)  into 
(2.49)  to  obtain  the  expression  for  the  disturbing  potential 
T  at  a  space  point  P  from  gravity  and  topography  data 
on  the  earth's  surface,  without  neglecting  the  topography  of 
the  earth  or  its  inclinations.  Only  the  second-order  deriva¬ 
tives  of  the  elevation  are  neglected  under  a  "smooth-earth” 
approximation . 

^  tco^sB 

y Ty^T^^^“‘^y 

O  X 

_ _ „  tan  $  taniv  ^  jo 

+  cosB  - ^ - * - }  dS  (2.51) 


If  we  set  the  inclination  components  t„  and  ly 
in  (2.51)  equal  to  zero,  and  use  (2.28-c),  then  (2.51) 
yields  the  following  spherical  formula  (see  Heiskanen  and 
Moritz,  1967,  p.l2): 

'’^P  •  TT  4n  ^  (2.5 


where  6g  is  now  the  gravity  disturbance  on  the  surface 
of  the  sphere  (S')  .  In  this  case,  (2.51')  corresponds 
to  the  "classical"  approach,  where  the  topography  is  not 
considered. 


It  now  remains  to  differentiate  (2.51)  with  respect 
to  X,  Y,  Z,  and  then  to  rotate  these  components  to  the 
(X,Y,Z)  system.  This  will  be  done  in  the  next  section. 


2.5  Differentiation  of  the  Disturbing  Potential 

In  order  to  obtain  the  derivatives  , 

we  first  compute  the  partials  — I  »  in  the 

(X,  7,  Z)  system  (defined  in  3X  SY 
section  2.2).  From  2.51  we  obtain: 


3  T 
3Xp 


3  T 
3  Zp 


1  T  tant  Jr  ^  (±)  dS 

4  ^  ■  cos6  r  X  3  Xp 


"  4^  jsrs  il;  <5^  <t>> 

■  t*'-.--  1  Ids 


3Xn  '3  y  i 


(2.52-a) 


=  ^  tani  l^(^)ds 

4  TT  ‘COs;i  r  x  Syp  < 


+  4-  (4:(4-))  tan-,  ])  dS 


’Yp 


3y'  i 


(2.52-b) 


=  i.  f  r  I  T 

4tt  •'gj  'cose 


cos  BtanQ 


tani. 


}  (4)  dS 


3Zp'  ^ 

^  T^y  ^  Vo'ieo^  ^^31^^  (^))tani^ 

^  tanty]  }  dS 


Zp 


(2.52-c) 


We  will  now  derive  one-by-lne  all  the  partial  deriva¬ 
tives  which  appear  in  the  equations  above.  The  spatial 
distance  il  between  the  space  point  P  and  any  surface 
point  Q  in  the  (X,  7,  Z)  system  is: 

A  =  ((Xp-Xq)2  + 


(2.53) 


From  (2.53)  we  obtain: 
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3 

Xp  -  Xq 

* 

3 

7p  -  ?o 
=  -  -^3^ 

3 

^  _  ^P- 

£' 

(2.54) 

We  also  have: 

— ^  ( 
3X 

-  -  siL 

2p"  3X 

— ^  ( 
ay 

=  .  1  ail 

51^  3y 

— ^  ( 
3Z 

=  -  1  Ml 

3Z 

(2.55) 

where 

3«.^ 

an^ 

aXp  3£^ 

“5x  5yJ 

3?p  ^  3£2  3Zp 

3x  agp  3x 

Ty  * 

an^ 

3Xp  ^  3£^ 

-w  s?;- 

3?p  ^  3£^  3Zp 

“sT  ^  ^ 

3Z^ 

3z 

3i^ 

55^ 

3Xp  ^  3£^ 

-if  *  ^ 

37p  .  dZ^  3^ 

~5z  51^  1z 

(2.56) 

or  in  matrix  form: 


Again 


from 

(2.53) 

for 

a  point 

3Xp 

=  2 

(Xp 

-  Xq) 

3Yp 

=  2 

(Yp 

-  V 

3il^ 

3  Zp 

=  2 

(2p 

- 

Q 
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fixed  on  S  we  obtain: 


(2.58) 


The  partial  derivatives  in  (2.57)  can  be  computed  from 
the_transformation  equation  between  the  (x,y,z),  and  the 
(X,Y,Z)  coordinate  systems: 


X  -  Xp 

X 

V  -  Yp 

=  Ra(-O0“-  ip))  R3(-AX)  R2(90“-  ig) 

y 

Z  -  Zp+rp 

(2.59) 

where  »  <tiQ  ^.re  the  geocentric  latitudes  of  the  points 
P  and  'Q  respectively,  and  rp  ,  rg  are  the  geocentric 
radii  to  these  points  (see  figure  2.17.  The  term  AX  is 
the  difference  in  longitude: 


AX 


(2.60) 


The  matrix: 


R  =  R2(-(90'’-  $p))  Rj(-AX)  RzOO®-  0g) 


in  equation  (2.59)  can  be  easily  shown  to 
sin^p  cosAXsin^p  -sinip  sinAX 

+  COS0p  COS^Q 


be  equal  to 

-sin^p  cosAX  cos^q 
+  cos$p  sin^Q 


R 


sinAX 


sin^Q 


cosAX 


-sinAX  cos$ 


Q 


-cos^p  cosAXsin^Q  cosip  sinAX 
+  sinig  cosiq 


cosip  cosAX  cosig 

+  sinip  sinig 


(2.61) 


3  /  ^  ^  \  "X 

37p  «- 

( _L  ( JL) ) 

^  ^3y  ^ 

-4-  (— 2-(^)) 

axp  32  «, 

-i-  (-2-(— )) 

avp  32  )i 

-i-  (-^(-i)) 

3Zp  32  2, 


=  ’  (%2e: 


■  3  Yp  X, 


■dy'  £ 


3  ,Xo-Xp^.„.  3  3  ,  Zo-Zp^ 


=  w  -2-  ( V^fi+w.  -2_(Ifiz|£)+w  ,-i-(-^a;^) 

^^3Xp^  27^^  ^^3Xp^  £3  ^^3Xp'  £^ 

(2.65) 


It  can  be  easily  shown  that  the  nine  derivatives  in  the 
right-hand  side  of  (2.65)  are  equal  to: 


(2.66) 


Let  us  now  substitute  the  partial  derivatives  (2.54), 
(2.65),  and  (2.66)  into  (2.52).  Since  we  want  to  evaluate 
thP  .'r.inponpn<  «  of  ^  «t  the  space  point  P  ,  let  us  also 
set  Xp  =  Xp  »  Zp  -  0  .  From  i2.52)  we  outain: 

+  3  ..  ^  W,,(5C^-  ^•)-kw3,X?  +  W3  3^Z 

^  S  cos  S 

-  2  cose  [(Wii(X2-  1^)  +  W,,X?  +  W13XZ)  tanXj^ 

+  (w,,  (3(2- 4-^)  +  WjzX?  +  W23XZ)  tanr  ]}  dS 
“  ^  y(2.67-a) 


1  =  ^  f  f  / — is 

3Yp  ^■’S 


tanT  } 


X  J 


3  ff  T  w,.YX  ■>•  w„  (?^>  "h  +  W3  3?2 

4rr  y  ^  S 

-2  cos  6[(Wii?X  +  W3  2(?^-  -|i)  +  WijYZ)  tantj. 


+  (Wji?X  +  W22(?2 - 1-)+  WjjYZ)  tanxy]  }  dS 

(2.67-b) 


^  .  1  r  f  {  _14  +  T  ££S|i^tanT  }  -|,  dS 

3Zp  ^  y  ‘  cose  r  X  '  *■  2 

3  f  f  T  r  W^iZX  H-  W32Z?  +  W33(Z^-  4^) 
■^57  y  ^  COSB 


-  2  coss  [(WiiZX  +  WizZ?  +  Wi3(Z2-  ^))  tanXj^ 

+  [W2iZX  +  W2  2Z?  +  W2  3(Z2-  -^  )  )  t3.nTy3}  dS 

(2.67-c) 


Finally,  rearranging  the  terras  in  (2.67)  we  arrive  at  the 
following  equations: 
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3T  _  1  5k  .  cos8  tan?  _ _  i  X 

JV  ^  - ? -  '“"'xl  T'  ‘>® 

+  j|-  //-p  (<i!*-  ^)  2cosB(WntanT^+w,itanTy)I 

s 

+  XY  [^;^  -  2cos3(Wi2tanTjj  +  W22'-anTy)] 

*  *  "'”‘‘”'y>"f2.68-a) 


3T  _  1  ,,,  +  r  cosBtan?  ,  ? 

T^  -  ^  V  ir- 


+  4I-  -  2cose(w,,tanT^  +  W22tanTy)] 

+  (?^-  -  2cos3(Wi2tanT^  +  WjztanTy)] 

■"  2':osB<»i.tanT,  * 


JT  .  1  rr(— T  tanTj  -Ir^s 

3Zp  47  Uoss  r  xJ  V 

*  ?l“  2cos6(WntanT^  +  W2itanTy)] 

+  Z1?  [gol^  -  2  cosg(Wi  2tanTjj  +  W22tanTy)J 

+  (Z^-  -|^)  [  -  2  cos6(wi  atanx^  +  W2  jtanTy )  ] }  dS 


3  '  '^cosB" 


(2.68-c) 


The  differential  surface  element  dS  in  all  the  equa¬ 
tions  above,  is  given  (in  general  form  for  an  ellipsoidal 
surface)  in  (Hotine,  1969,  eq.  30.72)  as: 


dS  *  (N  +  h)  (M  +  h)  cos())  d?  dX 


(2.69-a) 


where  N  ,  M  are  the  principal  radii  of  curvature  along 
the  prime  vertical  and  the  meridian,  and  h  is  the  height 
above  the  reference  ellipsoid.  Under  the  spherical  approx¬ 
imation  of  our  solution,  this  equation  becomes: 


dS 
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- 3-  cos  5  d  I  dX 

cos6 


( 2.69-b) 


If  the  surface  element  is  defined  by  the  intersections 
of  parallels  and  meridians  on  the  reference  sphere,  (2.69-b) 
reduced  to: 

“  cosT  ^^EAST  ~  ^WEST^  ^NORTH  “  ^SOUTH^ 

(2.69-c) 

The  data  that  is  needed  for  the  application  of  equations 
(2.68-a,b,c)  is  the  following.  (a)  Gravity  disturbances 
((Sg)  and  disturbing  potentials  (T)  on  the  surface  of  the 
earth,  (b)  the  elevations  of  the  surface  points,  and 
(c)  the  two  components  (tx  >  ty)  of  the  surface  inclinations. 
In  the  following  chapters  we  will  describe  how  this  data 
is  computed  for  our  simulation  studies. 


Rotation  to  the  Geocentric  System 


The  derivatives  (2.68)  are  referred  to  the  cartesian 
coordinate  system  (X,  Y,  Z)  located  at  the  space  point  P  . 
They  are  related  to  the  components  of  the  deflection  of 
the  vertical  and  to  the  gravity  disturbance  6g  at  P 
through  the  equations  (2.28-a,b,c) : 


3T 


-  yC 

+  Yh 

-  <5g 


(2.70-a) 
(2.70-b) 
( 2 . 70-c ) 


If  we  wish  to  transform  these  partial  derivatives  to  a 
geocentric  system  (X,  Y,  Z),  we  start  from  the  transformation 
between  (X,  ?,  2,),  and  (X,  Y,  Z): 

1  * 

?  -  Rz  (-(RO®  -  $p))  Rj  (-(180“-  Xp)) 

2 

( 


X 

Y 


Z 


or 


2  +  rp 

where 

-sinip 

cosXp 

-sinip  sinAp 

cosi. 

sinA 

P 

-cosAp 

0 

cosip 

cosAp 

cosip  SinAp 

sini. 

27. 

X 

Y 

Z 


(2.72) 


The  transformation  of  the  partials  in  (2.68)  to  the  geocen¬ 
tric  system  is  based  on  the  following  set  of  equations: 


3T 

3X  . 

3T 

3?  . 

3T 

3Z 

3X 

“TJT 

-IT 

3T  _ 

3T 

32  . 

3T 

3Y  . 

3T 

32 

5Y 

3X 

3Y 

“ST  ^ 

3T 

3T 

3T 

3?  . 

3T 

3Z 

3Z 

4. 

3Z 

3y 

3Z 

“ST 

or  in  matrix  form: 


3T 

3T 

-IT 

3T 

_  ,  ,(X,?,2), 

3T 

-  ^  (■nco)> 

3T 

3T 

-IT 

(2.73) 


where  the  Jacobian  matrix  J  is  computed  from  (2.71)  as: 


J 


,  (X,?,2)  , 


-sinip 

cosAp 

SinAp 

cosip 

cosA 

s 

-sinip 

SinAp 

-CosAp 

cosip 

sinX 

cos^p  0  sin0p 


(2.74) 
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In  summary,  equations  (2.68)  have  been  derived  for 
the  evaluation  of  the  derivatives  of  T  in  the  local 
(X,  Z)  system,  while  (2.73)  with  (2.74)  can  be  used 
to  compute  the  partials  of  T  in  the  geocentric  (X,  Y,  Z) 
system. 


2.7  Computation  of  the  Gravity  Vector. 

After  the  partial  derivatives  of  T  have  been  rotated 
to  the  geocentric  (X,Y,Z)  system,  and  the  gravity  distur¬ 
bance  vector  has  been  determined: 

5  -  grad  T  -  ,  -^)  ^  ^2. 3) 

the  components  of  the  gravity  vector  g  are  given  from 
equation  (2.5): 

n;  -  t  ^  /  3T  ST  ST  ,  ,  .  3T  3T  3T 

g  -  grad  V  -  I  +  y  - 

(2.3) 

U  being  the  normal  gravity  potential  of  the  equipotential 
of  the  equipotential  ellipsoid.  The  expressions  for  the 
partial  derivatives  of  U  with  respect  to  X,Y,Z  are  given 
in  (Heiskanen  and  Moritz,  1967,  chapter  6). 

For  some  purposes  we  need  the  vector  of  gravitation  f  , 
that  is  the  vector  of  the  pure  attraction  of  the  earth 
without  the  centrifugal  force  f  »  grad®  (Heiskanen  and 
Moritz,  1967,  p.228): 

f  =  grad  V  =  grad  (W  -  it)  =  g  -  grad*!*  =  g  -(uj^X,  aj^Y,0) 

(2.75) 


The  vector  g  contains  the  effect  of  the  rotation 
of  the  earth,  while  the  vector  f  does  not.  The  vector 
of  gravitation  is  of  considerable  interest  in  Geodesy, 
because  it  is  the  effect  of  the  gravitation  potential 
(V  =  W-<&)  of  the  earth  at  a  space  point,  which  is  not 
affected  by  its  rotation  (for  example,  an  artificial  satel¬ 
lite  which  is  orbiting  around  the  earth). 


Chapter  3 
THE  TERRAIN  MODEL 

3.1  Introduction 

In  order  to  avoid  the  errors  which  exist  in  real  data 
and  for  the  purpose  of  testing  equations  (2.68-a,b,c)  for 
the  gravity  disturbance  components,  we  decided  to  use  a 
simple  terrain  model,  with  synthetic  data  on  its  surface. 

As  it  will  be  described  later  in  this  chapter  two  disturbing 
masses  are  located  below  this  model  on  its  axis  of  symmetry, 
and  from  those,  the  exact  gravity  disturbance  vector  can 
be  computed  at  any  space  point.  The  gravity  disturbance 
6g  ,  and  the  disturbing  potential  T  can  also  be  computed 
at  any  point  on  the  surface  of  the  model.  This  data,  and 
the  inclinations  of  the  model  at  these  points  are  then 
used  in  (2.68)  to  compute  the  three  components  of  if  at 
any  space  point  P  ,  and  to  compare  them  with  the  exact 
components  evaluated  from  the  model's  disturbing  masses. 

The  idea  of  using  this  type  of  simulation  tests  is 
not  new,  and  the  details  will  follow  along  the  lines  of 
this  chapter.  For  gravity  anomaly  and  deflection  of  the 
vertical  computations,  a  conical  mountain  was  used  by 
Molodensky  et  al.,  (1962,  p.217),  andbyReit,  (1966). 

For  a  slightly  more  general  case  (the  model  is  on  a  sphere 
rather  than  on  a  plane  as  above),  see  Sjoberg  (1975,  p.82). 
Bjerhammar  (1963)  has  used  another  type  of  model  which 
is  a  homogeneous  spherical  cap  on  a  spherical  earth.  This 
model  has  also  been  discussed  by  (Sjoberg,  ibid,  p.79). 

In  order  to  find  estimates  for  the  effect  of  the  topographic 
elevations  on  the  external  gravity  anomalies,  Moritz  (1965) 
has  used  two  models  (a  conical  mountain,  and  a  two-dimensional 
model  e.g.  two  planes  intersecting  each  other).  Other 
simulations  were  performed  by  Koch  (1967-a,  1968-a)  with 
a  pyramid  and  a  cone.  These  models  are  similar  to  those 
mentioned  above,  but  they  have  the  advantage  of  allowing 
the  model's  inclination  to  vary  within  certain  limits. 

Among  the  simulations  above,  those  by  Koch  and  by 
Moritz  test  the  accuracy  of  the  computed  quantities  following 
the  Green's  identity  approach  The  rest  are  studies  on 
the  accuracy  in  the  computation  of  surface  or  space  quantities 
from  one  of  the  Bjerhamraar's  discrete  methods.  In  this 
paper  we  are  dealing  with  both  approaches.  The  simulations 
for  testing  the  Green's  approach  are  reported  in  chapter 
4,  and  those  for  the  discrete  (Dirac)  approach  are  reported 
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in  chapter  6.  In  the  present  chapter,  the  geometric  and 
the  dynamic  characteristics  of  the  model  are  described. 
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3.2  The  Geometry  of  the  Model 

The  terrain  model  which  has  been  selected  for  this 
study,  consists  of  a  cone  on  a  spherical  earth,  with  its 
top  replaced  by  a  spherical  cap  of  certain  radius,  such 
that  the  spherical  and  the  conical  surfaces  are  tangent 
to  each  other  at  their  intersection  (figure  3.1).  In  theory, 
this  rounding  of  the  model's  top  is  a  necessity  to  satisfy 
Lyapunov's  conditions  (Molodensky  et  al.,  1962,  p.217). 

Assume  that  the  axis  of  the  cone  intersects  the  surface 
of  the  mean  earth  sphere  (R  =  6371  km)  at  the  point  (0o  » 

A  0 t  h=0),  and  that  it  passes  through  the  center  of  mass 
0.  Let  (Oq  ,  »  h_)  be  the  coordinates  of  any 

point  on  the  surface  of'^the  model.  The  height  hq  is  the 
distance  from  Q  to  the  sphere  along  the  radius^from 
0  to  Q  .  Let  also  (ij>p  ,  Xp  ,  hp)  be  the  coordinates 
of  the  space  point  P  where  the  gravity  disturbance  vector 
is  to  be  computed.  At  the  point  Mj  and  on  the  axis, 
there  are  two  point  masses  mi  and  mj  respectively, 
which  generate  the  disturbing  potential.  Their  distances 
from  the  surface  of  the  sphere  (point  Ao )  are  ai  and 
a.2  •  As  inclination  of  the  model  (i),  we  define  the  angle 
between  its  surface  and  the  plane  tangent  to  the  sphere 
R  at  Ao  . 

As  we  will  show,  there  are  only  three  parameters  which  are 
needed  to  completely  describe  the  geometry  of  the  model 
in  terms  of  size  and  shape: 

1.  The  height  of  the  cone's  vertex  above  the  spherical 
earth  (h^) 

2.  The  extent  of  the  cone,  in  terms  of  the  angle  0  ,  and 

3.  The  radius  of  the  spherical  cap  (p). 

All  the  quantities  involved  with  the  model's  geometry 
are  computed  below  as  functions  of  the  three  parameters 
above.  Let  3  be  the  angle  at  the  geocenter  0  between 
the  axis  and  the  radius  to  point  D  where  the  cone  and 
the  spherical  cap  meet  each  other.  Let  also  u  be  the 
angle  at  the  geocenter  0  between  the  axis  and  the  radius 
to  an  arbitrary  point  Q  on  the  model.  There  are  three 


possible 

:  cases 

for  the  location 

of 

Q 

on 

the 

model : 

Q 

on 

the 

spherical  cap: 

0 

£ 

0) 

0) 

Q 

on 

the 

conical  surface: 

io 

0) 

< 

n 

Q 

on 

the 

mean  sphere  R  : 

U) 

180“ 

0)  can  be  computed  from: 


32. 


cosu)  ■  cos6o  COS0Q  3in9o  sine^  cos^X 

(3.1) 

where 

6o  *  90®  -  ^0 
-  90“  -  5,5 

AX  “ 

From  figure  3,1  we  easily  obtain: 

IB  .  hj  .  a  (1  -  oosn)  (3  2, 

and  the  inclination  i  of  the  model  is  then  computed  as: 


tani 


aB  _  ha  +  R(1  -  cosQ) 


BC 


R  sin  0 


(3.3) 


Then,  the  height  hy  of  the  model's  vertex  V  is  given 
by: 


(AK-p)=h^+p 


cosi 


(3.4) 


and  since 


5D  =  p  tani 


(3.5) 


we  obtain: 


tan  (i) 


p  si^ 

R  +  -  AD  sini 


(3.6) 


The  evaluation  of  the  height  of  any  surface  point 
is  done  according  to  its  location  on  the  model  (see  the 
three  conditions  above). 

Case-A  Point  Q  on  the  surface  of  the  cone  (w  i w i  n  , 
figure  3.1).  Applying  the  law  of  consines  to  the  triangle 
AOQ  we  obtain: 


^  +  i  -  0) ) 


which  gives: 
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Q  cosa)+  tani  sihw 

and  hence: 


(3.7) 


(3.8) 


The  distances  of  the  point  Q  from  the  two  masses 

are : 


vi  =  (Ri  +  Tq  -  2Ri  Tq  coso))^ 

and 

r2 

where 

Rx  =  R  +  ax 
R  2  “  H  “  21 2 


=  (R!  +  -  2R2 


’’Q 


COSU)  )  ' 


(3.9) 


(3.10) 


Case-B  Point  S  on  the  spherical  cap  (0*  <  w  <  ui  ,  figure 
3.1).  From  the  triangle  OKS  we  obtain: 

+  OK  -  2rg  oR  COSO)  ^2 

which  yields: 


rg  =  OE  COSU)  +  -  ijR*  sin^u) 

where 

OR  =  R  -t-  h^  -  p 

and  hence: 

hg  =  Tg  -  R 


(3.12) 


(3.13) 


(3.14) 


Case-C  Point  Q  on  the  surface  of  the  mean  sphere 
(n  <  u)  i  180®).  In  this  case 


34. 


h 


Q 


0 


(3.15) 


and  hence: 


r 


Q 


R 


(3.15-a) 


3.3  The  Data  on  the  Surface  of  the  Model 


After  the  geometry  of  the  model  has  been  established  in 
terms  of  its  size  and  shape,  the  disturbing  field  generated 
on  its  surface  and  in  the  exterior  space  is  now  determined. 

Four  more  parameters  are  needed  to  describe  the  dynamic 
characteristics  of  the  model.  These  parameters  are: 

1.  The  distances  ai  and  a2  of  the  two  masses  mi  and 
m2  from  the  surface  of  the  sphere  R  ,  and 

2.  The  gravity  disturbances  5gi  and  6g2  produced  by 
two  masses  at  the  vertex  V  of  the  model. 

Therefore,  for  the  simulation  studies  that  we  will 
conduct  in  this  paper,  we  need  to  define  the  following 
set  of  seven  parameters: 

h^  ,  ,  0  ,  ai  ,  a2  ,  6gi  .  ^gi 

The  disturbing  potential  at  any  point  Q  on  the  surface, 
or  in  the  exterior  space  is  computed  from: 

™  _  kmi  km  2 

Q  ~  ri  Tz  (3.16) 

where  ri  and  ri  are  given  by  (3.9),  and  kmi  ,  kma 
are  determined  from  the  given  boundary  conditions  6gi  , 

and  6g2  as  we  will  show  below. 

The  normal  component  of  the  gravity  disturbance  at 
any  surface  point  Q  is  computed  from; 

(3. IT) 

Using  (3.9),  (3.16),  and  (3.17)  we  obtain: 

_  kmi(rQ-R  cos^^)  .  km2(ro-R  cosw) 

r^ -  - ’■‘T] -  (3.18) 
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The  geocentric  distance  Tq  in  (3.18)  is  computed 

from: 

1.  equation  (3.7),  for  Q  on  the  conical  surface,  or 

2.  equation  (3.12),  for  Q  on  the  spherical  cap,  or 

3.  equation  (3.15-a),  for  Q  on  the  mean  sphere  R  . 

The  terms  kmi  and  km2  are  determined  from  the 
solution  of  (3.18)  using  the  boundary  conditions  6gi  and 
6g2  at  V  .  In  order  to  determine  luDi  and  kmz  ,  we  have 


kmi  =  6gi  (R  +  hy  -  Ri)^ 


and 


kmz  =  6gz  (R  +  h^  -  Rz)^ 


(3.19-a) 


(3.19-b) 


3.4  The  Gravity  Disturbance  Vector  in  Space.  Computed  from 
the  Model 

It  now  remains  to  find  the  expressions  for  the  compon¬ 
ents  of  the  gravity  disturbance  vector  as  a  space  point, 
which  is  generated  from  the  two  disturbing  masses: 


I  «  grad  T  =  ( 


3T 

3X 


3T 

3Y 


3T  , 

TZ' 


(3.20) 


The  disturbing  potential  T  is  given  by  (3.16),  where 
instead  of  (3.9)  we  now  use: 


ri  =  [(X 

-  Xi)"  + 

(Y  -  Yi )^  + 

(Z  -  Zi)'  ]* 

(3.21-a) 

and 

rz  =  [(X 

-  Xz)'  + 

(Y  -  Yz)'  + 

(Z  -  Zz)' 

(3.21-b) 

where  (X,  Y,  Z)  are  the  geocentric  cartesian  coordinates  of 
the  space  point  P  (5p  ,  Ap  ,  hp),  and  (Xi  ,  Yi  ,  Zi  ), 
(Xz  ,  Yz  ,  Zz),  are  the  cartesian  geocentric  coordinates  of 
the  two  disturbing  masses  mi  ,  and  mz  respectively: 

Xi  =  Ri  cos^o  cosAo 
Yi  *  Ri  cosio  sinAo 
Zi  *  Ri  sinjfl 


(3.22-a) 


aad 


36. 


X2  =  R2  cos$o  cosXo 
Y2  *  R2  cos^o  sinXo 
Z2  ■  R2  sin$o 


{3.22-b) 


Differentiating  (3.16),  and  taking  (3.21-a,b)  into 
account  we  obtain  the  following  expressions  for  the  three 
components  of  the  gravity  disturbance  vector  P  ,  as  com¬ 
puted  from  the  two  masses: 


,  dT  y  _  .  km  1  (X-Xi ) _  am2  ka.~a2  1 _  . 

P  [(X-Xi)2+  (Y-Yi)^+ (2-Zi)^ft  (X-X2)^+  (Y-Y2)?+  (Z-Z2)]^  = 


km2  (X-X2) 


(3.23-a) 


.  5T  )  =  km,  (Y-Y,) _  km2.(Y-Y^ -  ^ 

3y  ^  [(X-)4)^+ (Y-Yx)^+  (Z-Zx)^]'2  [(X-X2)^+  (Y-Y2)'+  (Z-Zz)^ 

(3.23-b) 


3T  X  a  .  kmi  (Z-2i) _ kmz  (Z-Z2) _ 

"^2:  P  "[(X-Xi)^+ (Y-Yi)'+  (Z-Zi)?]'^*[(X-X2>^+ (Y-Y2/+ (Z-Z2)'r' 

(3.23-c) 

Equations  (3.23)  yield  the  three  components  directly 
from  the  model,  and  therefore  these  results  are  errorless. 

In  the  next  chapter  we  will  discuss  the  simulations  tests 
we  made  in  order  to  investigate  how  well  equations  (2.68-a, 
b,c)  agree  with  the  exact  equations  (3.23)  in  the  case  of  the 
model,  without  integrating  over  the  whole  surface  of  the 
sphere . 


3.5  The  Effect  of  the  Model's  Symmetry  and  Center  of  Mass 

on  the  Simulations 

The  terrain  model  described  above  is  rotationally  sym¬ 
metric,  and  the  two  point  masses  located  on  its  axis  obviously 
generate  on  its  surface  a  disturbing  field  (6g  and  T  ), 
which  is  symmetric  too.  For  the  simulations  which  are 
described  in  the  following  chapters,  four  such  models  are 
used,  and  the  seven  parameters  for  each  one  of  them  are 
given  in  table  3.1  (these  parameters  were  defined  in  sections 
3.2  and  3.3). 


Table  3.1 
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Parameters  Defining  the  Models  for  the  Simulations 


i=10®.543 

i=20®.040 

i=39® .976 

Spherical 

inclination 

inclination 

inclination 

Cap 

Model 

Model 

Model 

Model 

4100  m 

4100  m 

4100  m 

23682973.2  m 

Q 

12’ 

6.08’ 

2.64* 

45® 

P 

200  m 

200  m 

200  m 

49-^5838.5  m 

ai 

2000  m 

2000  m 

2000  m 

2000  m 

a2 

4000  m 

4000  m 

4000  m 

4000  m 

6g, 

50  mgal 

50  mgal 

50  mgal 

50  mgal 

<5g2 

100  mgal 

100  mgal 

100  mgal 

100  mgal 

Note:  The  inclinations  are  computed  from  equation  (3.3). 

The  three  conical  models  will  be  referred  to  as  10"-,  20®-, 
and  40®-  inclination  models. 

Note  that  the  large  values  for  h^  and  p  in  the  list 
of  parameters  for  the  spherical  cap  model,  were  selected 
such  that  equation  (3.4)  will  yield  a  height  hy  of  the 
model's  vertes  V  above  the  mean  sphere  R  (figure  2.1) 
equal  to  4099.8  meters,  i.e.  of  the  same  magnitude  as  the 
heights  of  the  conical  models.  The  last  four  parameters 
are  the  same  for  all  models,  and  hence,  there  is  a  total 
gravity  disturbance  of  150  mgal  at  the  vertices  of  the  models. 
In  figure  3.2  through  3.5  the  gravity  disturbances  and  the 
disturbing  potentials  are  plotted  against  the  spherical  dis¬ 
tance  of  the  surface  points  from  the  model's  axis.  The  el¬ 
evation  of  the  surface  points  above  the  mean-earth  sphere 
R  are  also  given  at  2-kra  intervals. 

In  the  derivation  of  equations  (2.68)  following  the 
Green's  approach  (chapter  2),  no  assumption  was  made  about 
the  location  of  the  earth's  center  of  mass  with  respect  to 
the  origin  of  the  (X,Y,Z)  coordinate  system.  The  origin 
0  of  this  system  was  placed  at  the  earth's  geocenter  (see 
figure  2.1),  but  there  is  no  theoretical  requirement  for 
such  a  choice.  The  transformat ioi  equation  (2.59)  from  the 
(x,y,z)  system  to  the  (X,7,Z)  system  has  been  derived  by 
using  a  "geocentric"  (X,Y,Z)  system  whose  origin  does  not 
have  to  coincide  with  the  true  center  of  mass  of  the  earth. 

In  other  words,  there  will  be  no  effect  of  the  model's  center 
of  mass  on  our  simulation  tests  with  the  Green's  approach. 

As  long  as  we  have  the  boundary  values  6g  and  T  on  the 
surface  of  the  earth,  and  the  elevations  /  inclinations  with 
respect  to  a  reference  surface  (the  surface  of  the  mean- 
earth  sphere  in  our  case),  equation  (2^68)  can  be  used  to 
evaluate  the  components  9T/3X  ,  9T/9Y  ,  9T/92  in  the 

topocentric  (X,?,2)  system  centered  at  the  space  point  P  . 


The  rotation  from  this  system  to  the  "geocentric"  (X^Y ,Z) 
system  may  be  done  in  order  to  express  all  vectors  6  in 
a  unique  system. 

It  must  be  pointed-out  that  due  to  the  symmetric  char¬ 
acteristics  of  the  models,  there  are  only  two  distinct  compon¬ 
ents  of  the  gravity  disturbance  vectors  t  =  gradT  ,  namely 
3T/3r  ,  and  (1/r)  3T/3v  ,  where  is  the  angle  between 

the  model's  axis  and  the  radius  to  the  space  point  P  where 
0  is  computed.  On  the  other  hand,  we  decided  to  derive 
the  equations  for  the  computation  of  ^  in  such  a  way  that 
theycan  be  used  without  modification,  or  re-derivation  for 
a  future  application  with  real-data  (when  this  kind  of  infor¬ 
mation  becomes  available).  This  is  why  we  have  expressed 
the  components  of  ?  (equation  2.68)  not  in  a  polar  (r,  ) 

system,  but  in  a  cartesian  one.  Nevertheless,  the  transfor¬ 
mation  of  these  components  from  the  cartesian  to  the  polar 
coordinate  systems  can  be  easily  done,  using  (2.73),  and 
then ; 


3T 

3T 

3r 

Tx 

_al 

(r,|,X) 

3T 

3<t> 

3Y 

-11. 

3X 

3Z 

Finally,  the  coragonent  (1/r)  dT/dip  is  computed  from  the 
components  3T/3<))  ,  and  3T/3X  as: 


r  r 


/  ^2  1  i 


(3.25) 


The  direction  of  this  component  is  defined  from  the  space 
point  P  to  the  point  on  the  model's  axis  where  the  plane 
perpendicular  to  the  radius  r  at  P  interjects  the  axis. 

_The  Jacobian  of  the  transformation  from  (X,Y,Z)  to 
(r,  i,A)  can  be  computed  from  the  transformation  equations 
between  these  two  systems: 


X  *  r  cos  ^  cos  X 
7  »  r  cos  5  sin  X 
2  *  r  sin  ^ 


(3.26) 


and  then: 


4->  OS 


J 


((X.Y.Z)  ) 
(r,^,A) 


cos$  cosA 
-r  sin$  cosA 
-r  cos$  sinA 


cos^  sinA 
-r  sin|  SinA 
r  cos$  cosA 


sin$ 
r  cos$ 

0 

(3.27) 


For  the  simulation  tests  which  are  described  in  the 
following  chapters,  the  components  of  "s  will  be  compared 
to  the  exact  components  (from  the  model)  in  the  polar  (r,  if;) 
coordinate  system,  in  order  to  have  results  independent  of 
he  location  (^o  ,  Ao)  of  the  model's  center  on  the  sphere 
.  In  summary,  these  components  will  be  transformed  from 
the  topocentric  (X,?,Z)  system  to  the  polar  system  (r,  ip) 
using  equations  (2.73),  (3.24),  and  (3.25). 


41. 


(Sg  T 
(mgal )  (rf/sec^) 


Figure  3.4;  Simulation  Data  on  the  Surface  of  a  40“- 
Inclination  Conical  Model 
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Figure  3.5;  Simulation  Data  on  the  Surface  of  a  Spherical 
Cap  Model 


i 


Chapter  4 

SIMULATION  TESTS  WITH  THE  GREEN'S  APPROACH 
4.1  Introduction 


The  idea  of  the  simulation  is  to  compare  the  exact  grav¬ 
ity  vectors  as  computed  from  the  model,  with  those  evaluated 
from  Green's  approach,  using  surface  data  also  computed  from 
the  model. 

There  is  a  number  of  factors  that  will  affect  this  agree¬ 
ment.  The  inclination  of  the  model's  surface,  the  altitude 
of  the  space  point,  the  "density"  of  the  observations,  or 
in  other  words  the  grid  interval,  the  extent  of  the  area 
of  integration,  and  the  geometry  of  the  model  itself  (size 
and  shape),  are  among  these  factors.  Throughout  this  chap¬ 
ter  we  will  discuss  the  simulations  that  were  done  to  inves¬ 
tigate  the  effect  of  these  factors  on  the  agreement  between 
the  exact  and  the  computed  vectors.  Each  simulation  is  a 
four-step  procedure,  stunmarized  as  follows: 

Step-1.  Compute  the  exact  components  of  the  gravity  dis¬ 
turbance  vector,  using  equations  (3.23),  (3.24), 
and  (3.25). 

Step-2.  Compute  the  normal  gravity  disturbance  component 

6g  and  the  disturbing  potential  T  on  the  surface 
of  the  model  and  everywhere  else  within  the  speci¬ 
fied  area  (at  the  center  of  the  blocks  which  are 
formed  from  the  specified  grid).  Equations  (3.18) 
and  (3.16)  are  used. 

Step-3.  Compute  the  components  of  the  gravity  disturbance 
vector  at  the  space  point  P  ,  in  the  (X,^,Z)  co¬ 
ordinate  system  from  ( 2 .68-a ,b ,c ) ,  and  the  data 
generated  from  step  2.  Then,  rotate  these  three 
components  to  the  geocentric  (X,Y,Z)  system  using 
(2.73).  Finally,  compute  the  radial  and  the  hori¬ 
zontal  components  of  o  in  the  polar  (r,  <p)  system 
using  (3.24),  and  (3.25). 

Step-4 .  Compare  the  gravity  disturbance  components  from 

steps  1  and  3  in  terms  of  their  percentage  relative 
difference: 
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%  rel.  error  in  3T/9r  = 


3T 

^5r  exact  3r  '^computed 


exact 


X  100 


1  3T 

r  3i;j 


|r  H  exact  r  3i|)  computed 


■  3iJ;  "Uxact 


X  100 


(4.1) 


Assume  that  the  cone-like  terrain  model  is  located  on 
the  surface  of  the  earth  (more  specifically,  on  the  surface 
of  the  mean  sphere  R  ),  with  its  axis  intersecting  R  at 
the  point  (<>j  ,  ).  The  parameters  which  define  the  three 

conical  models  are  listed  in  table  3.1. 


In  order  to  apply  (2.68-a,b,c)  in  the  case  of  the  conical 
model,  we  divide  the  sphere  below  the  model  into  blocks  which 
are  bounded  by  meridians  and  parallels.  The  grid  spacing 
determines  the  density  of  the  data.  The  radii  passing  through 
the  intersections  of  the  grid,  form  on  the  surface  of  the 
model  another  grid(see  figure  4.1).  If  Q'  is  the  center 
of  the  block  on  the  sphere,  the  radius  through  Q'  intersects 
the  surface  at  the  point  Q  .  Equations  (3.16)  and  (3.18) 
yield  the  disturbing  potential  T  and  the  gravity  disturb¬ 
ance  6g  at  Q  .  Depending  on  the  location  of  the  block 
on  the  model  (on  the  spherical  cap,  on  the  conical  surface, 
or  on  the  sphere),  equations  (3.8),  (3.14),  or  (3.15)  are 
used  to  compute  the  height  of  Q  above  the  mean  sphere. 

The  inclinations  and  ly  in  the  north-south,  and  west- 

east  directions  at  '  Q  ,  are  computed  using  the  height  of 
Q  and  the  heights  of  the  neighboring  blocks  to  the  north, 
and  to  the  west  of  Q  ,  according  to  equations  below: 


tani^  = 

X  ®  __ 


tanx 

y 


hwest  -  h 

s 

y 


(4.2) 


where  h 

^north,  ^west 


the  height  at  the  center  (Q)  of  the  block, 
where  the  inclinations  are  to  be  determined, 
the  heights  at  the  centers  of  the  blocks 
to  the  north,  and  to  the  west  of  Q  . 
the  distances  from  the  centers  of  the  blocks 
to  the  north  and  to  the  west  of  Q  ,  from  it. 


The  disturbing  potential  T  ,  the  gravity  disturbance 
6g  ,  and  the  two  components  of  the  inclination  (t  ,  t  )  are 
the  data  needed  in  equations  (2.68),  along  with  tne  coordinates 
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of  the  surface  points  Q  .  These  equations  yield  the 
three  components  of  the  vector  at  the  space  point  P  in  the 
(X,  Y,  2)  system  located  at  P  ,  and  (2.73)  rotates  them 
to  the  geocentric  (X,Y,Z)  system.  Then,  (3.24)  and  (3.25) 
yield  the  two  distinct  components  of  ^  in  the  polar  (r,  ip) 
system.  The  exact  and  the  computed  vectors  are  finally  com¬ 
pared  in  terms  of  their  percentage  differences  from  equations 
(4.1).  For  the  investigation  of  the  effect  of  the  height 
(hp)  of  the  space  point  P  ,  these  vectors  are  computed  at 
different  altitudes,  at  equally-spaced  points  above  the 
sphere,  from  the  model's  axis  to  the  east  of  the  model. 


4.2  Simulation  Tests 

A  number  of  tests  were  made  to  investigate  the  influence 
of  the  following  parameters  on  the  accuracy  of  the  computed 
gravity  vector  components: 

a.  The  inclination  of  the  model. 

b.  The  geometry  of  the  model  (cone  vs.  sphere). 

c.  The  grid  interval. 

d.  The  size  of  the  integration  area,  which,  along  with  the 
grid  interval  determines  the  "density"  of  the  data.  Note 
that  the  center  of  the  integration  area  coincides  with  the 
center  of  the  model. 

e.  The  altitude  of  the  space  point. 

For  all  of  these  tests  we  assumed  that  the  model's  axis 
is  located  at 

^0  =  45“  ,  Xo  =  250“ 


This  choice  is  arbitrary,  but  since  our  software  was  designed 
for  real  as  well  as  for  simulated  data,  we  had  to  select 
a  "site"  for  the  model  in  terms  of  latitude  and  longitude. 

In  table  4.1  the  specifications  of  15  tests  are  summarized 
for  easy  comparison.  The  space  points  are  10'  apart  from 
each  other,  at  the  same  altitude  for  each  test,  from  250“ 
to  251“  east  longitude.  Therefore,  there  are  seven  points 
where  the  components  of  t  are  computed  and  compared  for 
each  test. 

The  first  five  tests  are  referred  to  a  conical  model 
whose  inclination  is  10“ .  From  tables  4.2  through  4.4  we 
see  that  by  increasing  the  integration  area  from  2°x  2“  to 
3“x3“,and  then  to  8“x8“,the  errors  at  the  distant  space 
points  (away  from  the  model)  are  reduced  dramatically,  but 
the  errors  above  the  model  are  almost  the  same.  Larger  inte¬ 
gration  areas  do  not  improve  the  results  above  the  model 
because  of  the  local  characteristic  of  the  disturbing  potential 
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of  the  two  masses  at  its  center.  However  for  the 
distant  points,  larger  integration  area  is  required  such 
that  there  is  enough  gravity  data  below  these  points. 

By  increasing  the  altitude  of  the  space  points  from 
10  km  to  20  km,  the  errors  are  reduced  both  o’er  the  model, 
and  at  the  distant  points  (table  4.2  vs  4.5).  From  tables 
4.6  and  4.7  we  see  that  the  errors  for  a  10'’-inclination 
are  slightly  smaller  over  the  cone  than  the  errors  for  a 
20°-inclination.  Tests  6  through  9  are  referred  to  a  20°- 
inclination  model. 

All  the  results  indicate  clearly  that  the  errors  above 
the  model  are  larger  than  the  errors  at  the  distant  points. 
Originally,  we  thought  that  the  reason  might  be  that  we 
should  have  used  a  finer  grid  over  the  model,  than  over 
the  neighboring  area  on  the  sphere,  because  of  the  fact 
that  most  of  the  disturbing  field  is  close  to  the  model's 
masses.  However,  by  using  0.5'  grid  in  the  inner  0‘’.3x0°.3 
area,  rather  than  a  I'-grid  (table  4.7,  vs.  4.8)  through¬ 
out  the  entire  4°x4°area,  we  did  not  see  any  improvement 
in  the  results.  The  errors  were  not  significantly  reduced 
neither  with  more  fine  0.1 '-grid  in  the  inner  zone  (table 
4.7,  vs.  4.9). 

For  an  altitude  of  100  km  above  the  mean  sphere  (table 
4.10),  the  errors  in  both  components  are  smaller  than  55( , 
that  is  much  smaller  than  the  corresponding  errors  at  an 
altitude  10  km  (table  4.7). 

On  table  4.11  we  see  the  errors  for  a  40®-inclination 
model.  If  we  compare  these  errors  with  those  on  table  4.7 
( 20'’-inclination) ,  we  see  that  they  are  much  smaller.  This 
is  probably  due  to  the  fact  that  a  40°-model  has  a  more 
"local"  effect  on  the  gravity  vector,  than  a  20°-raodel. 

The  inclinations  Tx  and  t  for  the  40®-model  are  larger 
as  compared  to  those  on  a  20^-model,  but  for  a  smaller  geo¬ 
graphical  extent. 

In  order  to  further  investigate  the  errors  from  a  very 
smooth-topography-raodel ,  we  substituted  the  conical  model 
by  the  spherical  model,  whose  parameters  are  listed  in  table 
3.1.  These  parameters  imply  a  height  of  the  model's  vertex 
V  above  the  mean-earth  sphere  R  ; 

hy  =  4099.8  meters 

which  is  of  the  same  magnitude  as  the  height  of  the  vertices 
of  the  conical  models.  In  other  words,  the  topography  is 
now  a  sphere  of  radius  o  ,  at  4.1  km  above  R  .  From  table 
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4.12  we  can  see  that  the  errors  are  now  very  small.  The 
components  of  the  gravity  disturbance  vector  are  computed 
from  the  surface  data  on  the  sphere  p  ,  with  an  error  which 
is  smaller  than  1%.  These  results  clearly  indicate  that 
the  Green's  approach  works  very  well  for  a  smooth  topography. 

From  the  eleven  tests  just  described,  we  can  draw  the 
following  general  conclusions: 

1.  The  higher  the  space  point,  the  smaller  the  relative 
errors  of  the  components  of  'J  become. 

2.  The  errors  are  larger  over  the  model,  then  over  the 
distant  zones. 

3.  For  a  given  grid  spacing,  and  for  a  fixed  altitude, 
the  increase  of  the  integration  area  does  not  improve 
the  results  over  the  model. 

4.  The  use  of  a  finer  grid  in  the  inner  zone  (O'. 5,  or 
O’.l  vs.  I'-grid),  does  not  improve  the  results  over 
the  model. 

5.  The  errors  are  very  small  for  a  smooth  topography  (the 
case  of  a  spherical  cap),  or  for  a  very  large  but  local 
topographic  feature  (the  case  of  a  40®-inclination  cone). 

Because  of  the  very  local  characteristics  of  the  models, 
there  are  no  truncation  errors  (i.e.  errors  caused  by  neg¬ 
lecting  the  information  outside  of  the  working  area).  In 
addition,  there  are  no  errors  in  the  surface  data,  since 
both  6g  and  T  are  rigorously  computed  from  the  model. 

The  only  factors  that  affect  the  results,  and  cause 
the  errors  in  the  computed  vectors,  are  the  following: 

1.  The  approximation  of  the  integral  equations  by  finite 
summations,  and  evaluations  of  the  kernels  at  the  center 
point  of  the  elementary  areas. 

2.  The  neglect  of  the  second-order  variations  of  the  topo¬ 
graphy  (see  equation  (2.43)). 

These  factors  are  responsible  for  the  larger  errors 
right  above  the  models,  as  compared  to  the  errors  far  away 
from  the  model,  or  to  the  errors  from  the  very  smooth  topo¬ 
graphy  of  the  spherical  cap  model..  However,  as  we  will 
see  in  chapter  8,  the  errors  from  the  classical  approach 
(the  direct  integration  method),  where  the  topography  is 
neglected,  are  larger  by  a  factor  of  2  than  the  errors  from 
the  Green's  approach,  indicating  the  significant  improvement 
of  the  later  method  in  the  computation  of  'J  . 

The  last  four  tests  in  table  4.1  are  applications  of  the 
Green's  approach  with  a  very  limited  amount  of  data  on  the 
surface  of  the  model.  These  tests  will  be  used  for  comparison 
purposes  with  the  Dirac  approach  in  chapter  6.  Only  576 
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blocks  have  been  formed  on  a  I’-grid  within  a  0®.4x0°.4 
area  on  a  10°-inclination  cone.  The  components  of  ^  are 
computed  at  5,  10,  20,  and  100  km  above  the  sphere  R  . 

Note  how  large  the  errors  are  at  a  5-km  altitude  (table 
4.13).  Also,  comparing  table  4.6  with  4.14,  we  see  that  at 
10  km  the  errors  are  of  the  same  order  of  magnitude,  indi¬ 
cating  once  more  that  the  effect  of  truncation  is  not  signi¬ 
ficant  (for  these  two  tests,  only  the  size  of  area  is  dif¬ 
ferent  ) . 


4.3  Related  Work,  and  Comparison  with  Koch's  Method 

Among  the  papers  which  have  been  published  on  the  compu¬ 
tation  of  certain  components  of  the  gravity  disturbance 
vector  using  the  Green's  third  identity,  the  most  elegant 
approach  in  our  opinion  is  the  one  by  Koch  (1968-a).  There¬ 
fore,  it  is  essential  to  comment  on  the  similarities  and 
on  the  differences  between  his  and  our  approaches,  and  to 
discuss  his  results  under  the  light  of  ours. 

Koch's  approach  is  a  simulation  study  with  a  cone¬ 
like  model  as  in  our  case,  but  the  solution  for  the  disturbing 
potential  and  its  partial  derivatives  is  obtained  through 
an  iterative  procedure.  In  addition,  the  computed  quantities 
are  surface  quantities,  and  not  at  space  points.  Furthermore, 
his  equations  have  been  derived  specifically  for  the  simulation 
studies  that  he  performed,  and  allow  for  the  maximum  inclina¬ 
tion  of  the  cone  (B)  only  (not  a  two-component  ,  ty 

description  of  the  topography).  Therefore,  his  equations 
cannot  be  used  for  a  real-world  application  without  additional 
work.  Also,  the  whole  approach  is  based  on  a  planar  approx¬ 
imation  . 

The  large  errors  that  were  found  by  Koch  in  the  case 
of  inclinations  larger  than  20®,  we  believe  that  are  due 
to  the  fact  that  the  term  T(PoHs  included  in  his  integral 
equation  (21),  while  the  integration  is  not  carried  out  over 
the  whole  surface  of  the  earth  (according  to  the  discussion 
below  our  equation  (2.9')).  In  addition,  the  fact  that  Koch's 
equations  were  derived  for  computations  of  the  derivatives 
of  T  on  the  surface  of  the  earth  only,  does  not  permit 
comparison  of  our  results  with  his. 
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Computed  Rigorously,  and  From  Integration  Procedures  Using  the 
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Chapter  5 


DISCRETE  APPROACHES  FOR  THE  COMPUTATION  OF  THE 
GRAVITY  DISTURBANCE  VECTOR  IN  SPACE 

5.1  Introduction 

Throughout  the  three  previous  chapters,  we  developed  and 
tested  the  equations  for  the  computation  of  the  components 
of  the  gravity  disturbance  vector  at  a  space  point.  The 
described  solution  is  based  on  the  application  of  Green's 
third  identity,  the  data  being  given  on  the  physical  surface 
of  the  earth.  Hence,  the  method  follows  the  modern  approach 
for  the  solution  of  the  boundary  value  problem  for  the  deter¬ 
mination  of  the  external  gravity  field  of  the  earth.  Both, 
the  coventional  solution  to  the  geodetic  BVP  (using  the 
Stokes'  integral),  and  the  modern  solution  to  the  geodetic 
BVP  (following  Molodensky),  require  the  data  be  given  as 
mean  values  over  small  surface  elements.  Furthermore,  the 
conventional  approach  requires  gravity  reduction  to  the 
geoid,  a  process  which  assumes  knowledge  of  the  earth's 
densities,  while  Molodensky's  solution  does  not  require 
any  reductions. 

In  practice,  gravity  measurements  are  taken  at  a  finite 
number  of  stations  only.  The  need  to  solve  the  BVP  in  space 
from  a  finite  number  of  gravity  observations,  let  Bjerhammar 
to  develop  some  methods  (documented  in  a  number  of  papers 
which  are  referenced  below),  in  general  called  discrete  ap¬ 
proaches  to  the  solution  of  the  BVP  of  physical  Geodesy. 

The  formal  statement  by  Bjerhammar  (1963)  reads  as 
follows: 

"A  finite  number  of  gravity  data  (Ag)  is  given  for 
a  non-spher ical  surface,  and  it  is  required  to  find 
such  a  solution  that  the  boundary  values  for  the  gravity 
data  are  satisfied  in  all  given  points". 

For  a  given  finite  set  of  gravity  observations  (or 
more  specifically  gravity  anomalies  which  will  be  defined 
more  precisely  in  the  next  section),  there  is  always  a 
f ictitious  field  of  gravity  anomalies  on  an  internal  sphere 
that  satisfies  the  given  boundary  values  on  the  surface. 

This  sphere  is  completely  imbedded  inside  the  earth,  and 
appears  in  the  literature  under  the  names  "Bjerhammar  sphere", 
or  "geosphere".  Bjerhammar 's  methods  are  applications  of 
the  Poisson's  integral  equation  and  Stokes'  formula  (Bjer¬ 
hammar,  1976,  1978;  Sjoberg,  1978).  Earlier  solutions 
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59. 

considered  the  gravity  anomalies  on  the  geosphere  as  mean 
anomalies,  but  later  solutions  have  been  formulated  on  the 
basis  that  the  reduced  anomalies  are  point  anomalies  on 
the  geosphere  with  no  block  size  associated  with  them 
(Bjerhammar,  1978,  p.220).  The  common  characteristic  of 
all  methods  is  that  a  finite  set  of  point  gravity  anomalies 
is  given  on  the  physical  surface  of  the  earth 

A  special  case  of  Bjerhammar 's  reflexive  prediction 
method  (Bjerhammar,  ibid),  is  the  Dirac  approach,  which 
has  been  compared  to  collocation  by  Sjoberg  (1978),  This 
approach  is  based  on  the  computation  of  the  fictitious  gravity 
anomalies  at  the  so-called  carrier  points,  located  at  the 
intersections  of  the  radii  to  the  observation  points,  with 
the  geosphere.  The  radius  (R3)  of  the  geosphere  was  found 
to  be  a  critical  parameter  in  this  type  of  solutions.  As 
(Sjoberg,  1978,  p.64)  reported,  the  ootimum  depth  of  the 
geosphere  from  the  mean  earth  sphere  R  ,  is  approximately 
half  the  mean  distance  between  surface  neighboring  obser¬ 
vations. 

In  the  sections  that  follow,  we  will  review  the  theory 
behind  these  discrete  approaches,  with  emphasis  on  two 
methods,  namely  the  mean-value  approach,  and  the  Dirac  ap¬ 
proach.  Then,  the  Dirac  approach  will  be  applied  to  simu¬ 
lated  data  from  the  model  described  in  chapter  3,  and  the 
results  (the  components  of  the  disturbance  vectors)  will 
be  compared  to  the  exact  vectors  from  the  model,  just  as 
we  did  with  Green’s  approach  in  chapter  4. 


5.2  Theoretical  Development 


Let  us  start  from  the  statement  on  the  Dirichlet's 
problem  or  the  first  BVP  of  potential  theory  (Heiskanen 
and  Moritz,  1967,  p.34): 

"Given  an  arbitrary  function  on  a  surface  S  ,  deter¬ 
mine  a  function  V  which  is  harmonic  either  inside 
or  outside  S  ,  and  which  assumes  on  S  the  values 
of  the  prescribed  function". 

Dirichlet's  problem  can  always  be  solved  if  S  is 
the  surfr.ce  of  a  sphere,  an  explicit  solution  being  given 
by  the  Poisson's  integral  equation  (Heiskanen  and  Moritz, 
1967,  p.35),  which  for  the  exterior  of  S  is  written  as: 


60. 


V  is  an  arbitrary  harmonic  function  on  S  and  in 
the  exterior  space,  (9',  X')  are  the  polar  coordinates  of 
the  variable  point  on  the  sphere  (R),  (r,  6,  X)  are  the 
coordinates  of  the  exterior  point,  and  I  is  the  distance 
from  (r,  9,  X)  to  (R,  9*,  X'): 

i  =  (r^  +  R^  -  2R  r  cosi;^)^ 


cosili  =  cos9  cos9 '  +  sin9  sin9'  cos(X'  -X) 


The  function  rhg  is  a  harmonic  function  (Heiskanen 
and  Moritz,  1967,  p.90),  and  hence,  Poisson’s  Integral  can 
be  applied  for  any  point  in  space  as: 


rpSgp 


//  da 

o 


(5.4) 


do  =  sine*  d9  •  dX ' 


(5.5) 


is  the  surface  element  on  a  unit  sphere.  For  the  geosphere 
(Rg),  equation  (5.4)  becomes: 


where 

dS  =  Rg  do  (5.7) 

S  being  the  area  of  the  geosphere  B  (S  =  4tt  Rg).  The 
quantities  which  appear  in  (5.6)  are  defined  as  follows 
(see  also  in  figure  5.1): 

Agp  :  The  gravity  anomaly  at  the  space  point  P 

Agp  *  Kp  ”  Yp 

where  P'  is  that  point  along  the  vertical 
through  P  for  which  Up,  *Wp  (cf.  Sjdberg, 

1978); 

Ag*  :  The  fictitious  gravity  anomalies  on  the  geosphere, 
which  generate  the  surface  Ag  ,  and  which  have 
to  be  determined; 

r  ;  The  geocentric  distance  to  P  . 

P 


61. 


Obviously,  equation  (5.6)  also  holds  for  any  point 
Q  on  the  physical  surface  of  the  earth,  because  the  geos¬ 
phere  has  been  defined  to  be  completely  imbedded  inside 
the  topographic  masses,  and  therefore,  the  requirement  of 
the  Poisson's  integral  (5.1)  is  met  (the  solution  is  given 
for  points  outside  of  the  surface  S  ,  the  surface  being 
the  sphere  Rg  in  our  case).  Equation  (5.6),  applied  for 
a  surface  point  Q  ,  becomes: 


-.Rr 

4TTrQ 


// 

B 


By  substituting: 


A_g* 

tT 


dS 


(5.8) 
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and 


D,  =  — -  .  /I  -  2  eos»  > 

in  equation  (5.8),  we  obtain: 


(5.9) 


(5.10) 


(5.11) 


The  precise  definition  of  the  surface  gravity  anomaly 
AgQ  is  similar  to  the  definition  of  the  space  anomaly 
Agp  (see  below  equation  (5.6)).  In  other  words,  AgQ  is 
the  difference  between  the  measured  gravity  gQ  at  the 
surface  point  Q  ,  and  the  normal  gravity  yq  computed 
at  that  point  Q'  along  the  ellipsoidal  normal  through 
Q  ,  for  which  Uq*  •  All  points  Q'  form  the  surface 

of  the  telluroid  (Heiskanen  and  Moritz,  1967,  p.292). 

The  central  idea  in  the  Bherhammer's  methods  is  to 
apply  an  integral  equation  like  (5.8)  or  (5.11)  to  the  given 
surface  data  AgQ  ,  and  to  compute  the  fictitious  Ag*  on 
the  geosphere,  in  such  a  way  that  all  AgQ  are  satisfied 
by  the  analytically  determined  Ag*  .  Furthermore,  the 
gravity  anomaly  at  a  space  point  P  ,  determined  from  the 
anomalies  on  the  geosphere,  must  be  identical  to  the  gravity 
anomaly  computed  from  the  measured  surface  AgQ  ,  provided 
that  there  are  no  errors  in  the  data,  and  that^the  integrations 
are  performed  over  the  whole  surface  of  the  earth.  This 
is  a  consequence  from  the  uniqueness  of  Stokes'  theorem 
(Heiskanen  and  Moritz,  1967,  p.l7). 
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Let  us  now  show  how  an  integral  equation  like  (5.11) 
can  be  applied  to  a  discrete  approach. 


5.2.1  The  Mean-Value  Approach 


If  the  surface  of  the  geosphere  is  divided  into  a  finite 
number  (N)  of  blocks  (AS.),  equation  (5.11)  can  be  rigor¬ 
ously  written  as:  ^ 


Q 


= 

47rR2 


to) 


N 

I 

i®l 


/; 

ASi 


Ag* 

"cr 


dS 


(5.12) 


or,  approximately: 


(5.13) 


where  Ag*  is  evaluated  at  a  certain  (but  unknown)  point 
inside  A^.  ,  usually  assumed  to  be  the  center  of  the  block, 
in  which  case  (5.13)  is  an  approximation.  The  conventional 
approach  to  convert  a  theoretical  continuous  integral  trans¬ 
formation  like  (5.11),  to  a  discrete  summation  transforma¬ 
tion  like  (5.13),  is  called  discretization  (see  for  example 
in  Robertson,  1978,  p.4-1  ).  Discretizations  are  of  great 
importance  in  geodesy,  since  we  have  to  deal  with  finite 
sets  of  data,  irregularly  distributed  on  the  surface  of  the 
earth.  The  substitution  of  integral  equations  by  finite 
summations  is  a  common  practice,  especially  for  computer- 
oriented  applications.  Obviously,  the  best  discretization 
process  is  the  one  which  yields  the  minimum  error  of  the 
computed  quantity  with  respect  to  the  true  value.  Robert¬ 
son  (ibid)  discusses  three  discretization  processes,  and 
suggests  a  method  for  evaluating  their  accuracies  by  comparing 
their  spectra. 


The  mean-value  approach,  requires  the  evaluation  of 
the  areas  AS^^^  for  each  block,  and  from  the  theoretical 
point  of  view^(5.13)  holds  as  an  approximation  since  the 
location  of  the  point  for  evaluation  Ag*  is  unknown.  The 
Dirac  approach  which  will  be  discussed  later  avoids  this 
problem  because  it  determines  Ag*  at  the  carrier  points 
whose  positions  are  pre-deterrain4d,  and  therefore  are  precisely 
known . 


If  the  subdivision  of  the  geospher  in  blocks  is  such 
that  their  n\xmber  (N)  is  equal  to  the  number  of  observations 
AgQ  ,  then,  a  direct  solution  to  (5.13)  for  the  fictitious 
Ag^  ,  could  be: 


r^LS. 
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Ag* 


(5.14) 


where  the  underlines  mean  matrixes,  or  vectors,  and  the 
elements  of  A  are: 


t?i  (1  -  t?) 


aS  i 

Tfr 


(5.15) 


However,  such  an  approach  would  be  practically  impossible 
if  a  large  amount  of  surface  data  is  available.  Earlier 
solutions  to  the  Bjerhammar's  problem  used  an  iterative 
approach  to  compute  Ag*  from  Ag  (see  for  example  Heiskanen 
and  Moritz,  (1967,  p.318),  or  Emrick,  (1973)).  In  the  next 
section  we  describe  the  Gauss-Seidel  iteration  technique, 
and  we  investigate  the  possibility  for  employing  an  acceler¬ 
ation  procedure. 


5.2.2  Gauss-Seidel  Iteration  Method 

For  a  linear  system  of  n  equations  with  n  unknowns; 

-  -  “  -  (5.16) 
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(5.17) 

the  Gauss-Seidel  iteration  method  to  solve  is  (Carnahan, 

et  al.,  1969,  p.299): 

J  1+i  J«1  (5.18) 


where  (k)  denotes  the  iteration  step.  In  this  iterative 
method,  the  newly -computed  components  of  the  solution  vector 
X  are  always  used  in  the  right-hand  side  of  (5.18)  as  soon 
as  they  are  obtained.  A  sufficient  condition  to  guarantee 
the  convergence  of  the  iteration  is: 


<  1  . 


1  i  i  i  N 


''5.19) 


This  cone ition  will  later  be  used  to  coranute  the  opti- 
miun  radius  of  the  geosphere  in  order  to  guarantee  the  conver¬ 
gence  of  the  solution  for  Ag*  . 


As  an  example,  let  us  apply  (5.18)  to  the  linear  system 
(5.13).  Assume  that  N  surface  gravity  anomalies  Ag  are 
given,  and  that  we  want  to  determine  the  N  unknowns  Ag* 
on  the  geosphere.  The  surface  Ag  are  point  anomalies, 
but  the  Ag*  on  the  geosphere  are  mean  anomalies,  associated 
with  a  certain  block  size  which  depends  on  the  density  of 
the  surface  data.  Let  Q  be  the  points  on  the  surface 
where  the  data  is  given,  and  Q'  be  the  center  of  the  block 
on  the  geosphere  (see  figure  5.2). 


With  this  notation,  the  linear  system  (5.13)  becomes: 


t^  (1-  th 


4. 


(5.20) 


with 
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t  *  trt  = 

.)  Q,1  tqj 

The  application  of  the  Gauss-Seidel  method  (5.18)  to  the 
system  (5.20)  yields: 


(k+1) 

- 

where  i  =  1 


.  .  .  ,  N 


^5.21) 


k  = 


1.  2, 


and 


Ag* 


(‘t=l) 

Q'j 


maximum  number  of 
^^max> 


iterations  allowed 


j  =  1,  2,  .  .  .  ,  N  as  an  initial 
approximation 


An  equation  analogous  to  (5.21)  but  somewhat  simpler  will 
be  derived  later  for  the  Dirac  approach.  The  above  iteration 
for  the  unknows  Ag*  can  be  stoped  either  if  we  exceed 
the  maximum  number  of  iterations  allowed,  (say  kj^3^jj  =  30), 
or  if  two  successive  Ag*  at  all  points  (iteration  k  and 
k  +  1)  are  different  by  no  more  than  c  ,  a  quantity  which 
is  as  small  as  we  desire. 


Acceleration  Techniques.  Among  the  existing  methods 
for  the  iterative  solution  of  large  systems  of  equations, 
the  Gauss-Seidel  (’’successive  iterations")  method  converges 
about  two  times  faster  than  the  Jacobi's  ("simulation  iter¬ 
ations")  method  (cf.  Westlake,  1968,  p.56).  In  addition, 
there  ar  certain  acceration  techniques,  which  yield  a  much 
larger  rate  of  convergence  than  the  two  methods  mentioned 
above.  An  acceleration  technique  requires  the  splitting 
of  the  matrix  B  according  to  (Isaacson  and  Keller,  1966, 
pp. 73-80) : 


where  B  =  N(a)-P(a) 

N  (a)  =  (1  +  a)  No 
P  (a)  =  (1  +  o)  No  -  B 

and 
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The  estimation  of  the  optimum  parameter  Ca=cig  ^)  which 
yields  the  largest  rate  of  convergence,  requires  the  knowledge 
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of  the  minimum  and  of  the  maximum  eigenvalues  of  the  matrix 
Pq  ,  according  to  the  formula: 


opt 


IffllTL 


^max 


After  the  parameter  otQpt  has  been  determined,  the 
iterative  solution  to  the  system  ^=U  can  be  obtained 
as  (Westlake,  1968,  p.62,  where  u)=  1  +  a): 

.(k+1)^  1  +  g  r..  .  ?  b.  . 


^i 


^  j  =  i+l  J=1 


J  =  i-^i  ■  (5.18-b) 

which  for  a=0  reduces  to  the  Gauss-Seidel  technique  (5.18), 


ment 
will 
(a) . 


(b) 


Clearly,  from  an  application  point  of  view,  the  employ- 
of  an  acceleration  technique  requires  that  the  user 
follow  one  of  the  following  two  procedures. 

Form  the  matrix  N£^Pi  >  find  its  eigenvalues  X^jin 
and  Xjnax  »  then  determine  Oopt  which  yields 

fastest  rate  of  convergence  in  (5,lS-b). 

Emoirically  estimate  a  parameter  a  which  yields 
a  rate  of  convergence  faster  than  that  of  the  Gauss- 
Seidel  technique.  Such  an  empirical  determination 
(used  for  the  downward  continuation  of  gravity  anom¬ 
alies),  has  been  used  by  Koch  (1968-b).  Westlake 
(1968,  pp. 62-63)  describes  a  procedures  for  the  estim¬ 
ation  of  a  but  in  the  special  case  that  B  is  a 
symmetric  matrix.  Finally,  another  empirical  procedure 
for  the  estimation  of  a  is  described  in  (Isaacson 
and  Keller,  1966,  pp. 78-80),  but  the  method  is  graphical, 
and  hence  not  suitable  for  computer-oriented  applications. 

For  large  matrices,  the  analytical  procedure  (a)  above 
requires  considerably  additional  work,  before  the  actual 
iteration  (5.18-b)  starts.  Furthermore,  the  empirical  tech¬ 
niques  (b)  mentioned  above  are  either  not  fully  automatted, 
or  not  applicable  for  a  general  (not  symmetric)  matrix. 
Therefore,  and  for  reasons  of  simplicity,  we  decided  to 
use  the  Gauss-Seidel  iterative  method. 


5.2.3  The  Computation  of  the  Gravity  Disturbance  Vectors 
in  Space  (Mean-Value  Approach) 

Once  Ag*  have  been  determined,  the  evaluation  of 
the  three  components  of  the  gravity  disturbance  vector  in 
space  can  be  done  following  standard  procedures.  The  three 
components  of  T  >  are  defined  as  (Heiskanen  and  Moritz, 
1967 ,  equ.6-29) : 
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(Sgr^p 


(JJT) 

3  r'^P 


(5g;^)p 


1  /  3  T 

Tq  cos^_  “ST*  p 


(5.22) 


where  (rp  ,  $p  ,  Xp)  are  the  coordinates  of  the  space  point 
P  where  these  components  are  to  be  computed.  The  disturbing 
potential  at  P  is  given  by  Pizzetti's  generalization  of 
Stokes’  formula  (Heiskanen  and  Moritz,  1967,  equ.8-88): 


Tp  =  T(rp,  ^p,  X)  =  ^  J|Ag*  S(r,^:  do 

B 

From  (5«22)  and  (5.23)  follows  that: 

/jAg*  -  cosa  do 


(5.23) 


(5.24) 


Since  Ag*  represents  the  value  of  the  fictitious 
anomaly  over  a  block  on  the  geosphere,  we  apply  the  same 
approximation  as  we  did  in  section  5.2.1,  to  obtain  (Note: 
dS  =  R|da ) : 


‘Si 

”*x'p  '‘Si  ' 
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The  index  i  refers  to  the  center  of  the  block  AS.- 
this  is  in  approximation,  since  the  point  to  evaluate  the 
kernels  in  (5.24)  is  unknown  -  ,  and  a.  is  the  aximuth 
from  the  projection  of  P  on  the  geospfter  to  the  point  Q!  : 


tana^  = 


_ cosji  sin<Ai~Ap)  _ 

cos^p  sin“^  -  sin^p  cos$j^  cos(A^-Ap) 


(5.26) 


The  partial  derivatives  of  the  Stokes'  function 
respect  to  r  and  o  are  computed  as: 

cosi/'  =  cosi^/p^  =  sin5p  sin^^^^  +  cos^p  cos^^  cos(A^ 
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the  notation  above  we  have: 
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Finally,  the  components  of  the  gravity  disturbance 
vector  in  a  geocentric  coordinate  system  are  (Heiskanen 
and  Moritz,  1967,  equ.6-18): 
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5.3  The  Dirac  Approach 

The  Dirac  approach  avoids  the  continuation  of  the  sur¬ 
face  anomalies  to  mean  (fictitious)  anomalies  on  the  geos¬ 
phere.  Instead,  the  gravity  information  of  the  N  surface 
point  gravity  anomalies  Ag  ,  is  downward-continued  to 
gravity  anomaly  impulses  (spikes)  Ag®  at  a  finite  number 
of  carrier  points  (defined  in  section  5.1).  This  downward 
continuation  is  an  analytic  continuation,  by  which  is  meant 
that  the  N  spikes  satisfy  the  surface  data,  without  having 
any  physical  meaning  (i.e.  this  is  not  a  gravity  reduction 
that  yields  the  true  anomaly  on  the  geosphere). 

If  the  number  of  the  gravity  impulses  Ag®  on  the 
geosphere  is  equal  to  the  number  of  the  given  surface  gravity 
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anomalies  Ag  ,  the  case  is  called  "non-singular"  (Bjerham- 
mar,  1977),  and  the  carrier  points  can  be  selected  at  the 
projections  of  the  surface  points,  on  the  geosphere,  along 
the  geocentric  radii  (see  figure  5.3).  Qj^  are  the  surface 
points  where  the  data  is  given,  and  are  their  projections 
on  the  geosphere,  where  the  spikes  Ags  are  to  be  computed. 


Figure  5.3;  The  Dirac  Approach:  The  Spikes  Ag®  at  the 
Carrier  Points 


Let  us  now  define  the  anomaly 
(5.11)  as: 

N 

Ag*  =  I  Ag®  6(r  -  r.  ) 
k=l  ^ 


Ag*  which  appears  in 


(5.33) 


where  r  :  the  radius  to  the  current  point  on  the  geosphere. 

rj.:  the  radius  to  the  k^"  spike  Agg  on  the  geosphere 

and  Afr-rj,);  the  Dirac  "delta"  function  defined  through 
the  following  integral  equation  for  an  arbitrary 
function  f(r)  on  the  geosphere. 
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Substituting  (5.33)  into  (5.11)  we  obtain: 
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and  using  (5.34)  with  f(r)  =  1/Dq  i  we  obtain: 
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iSo  =  to  <l-‘o)  ^ 


(5.35) 


The  reader  can  see  that  equation  (5.35)  of  the  Dirac  ap¬ 
proach,  corresponds  to  equation  (5.13)  of  the  mean-value 
approach,  the  difference  being  that  the  former  does  not 
involve  any  areas  AS^  associated  with  the  spikes.  Further¬ 
more,  the  positions  of  the  carrier  points  are  precisely 
known,  and  hence  the  Dirac  approach  avoids  the  approximation 
of  the  mean-value  approach,  where  the  kernels  are  evaluated 
at  the  centers  of  the  blocks. 


5.3.1  Iterative  solution  for  the  Spikes  AgS 

The  computation  of  the  N  spikes  AgS  from  the  given 
finite  set  of  observations  Ag  ,  is  done  using  the  Gauss- 
Seidel  iterative  method  described  in  section  5.2.2.  More 
specifically,  from  (5.18)  and  (5.35)  we  obtain: 


equation  (5.36)  yields: 
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Ji^i 


"C- 


(5.37) 


As  initial  values  we  used  -  0,  j  =1,  2,  .  .  .  N, 

and  not  Ag^'.-  Ag^.  ,  because  fr^om  our  preliminary  tests 
we  found  that^  the  magnitude  of  the  spikes  is  much  smaller 
(or  the  order  of  10~^  ®m.sec~^ )  than  the  mamitude  of  the 
surface  anomalies  (from  10~^  to  10“®m. sec~^) . 


The  iteration  scheme  (5.37)  can  be  terminated  either 
if  it  exceeds  a  specified  maximum  or  if  two  successive 

iterations  for  a  spike  yield  values 


and 


(k+1) 


which  are  different  by  less  than  a  specified  limit.  For 
example,  if  we  want  to  have  a  12-digit  agreement  between 
these  two  values,  we  stop  the  iteration  whenever  the  fol¬ 
lowing  condition  is  satisfied  for  every  spike: 

.(k+1)! 


Ag' 


^.(k+1) 


,(k) 


-  Ag' 


Agf 


10 


1  2 


(5.38) 


5.3.2 


The  Convergence  of  the  Iterative  Solution,  and  the 
Radius  of  the  Geosphere 


Let  us  now  examine  the  conditions  for  the  convergence 
of  the  Gauss-Seidel  iterative  method  for  the  spikes.  A 
sufficient  condition  for  the  convergence  of  this  method 
was  given  in  section  5.2.2  as  inequality  (5.19).  By  applying 
this  condition  in  the  case  of  the  Dirac  spikes  (5.35)  we 
obtain : 


j  =  1,  2, 


which  becomes: 


N: 


N 

I 

i=l 


fi 


(1  -  t^  ) 

unit 


( 1  -  ) 

rrrt-t- 


<  1 


j  *  1,  2, 


N.  ^  ( 1  -  t.1  )  ^ _ ^  ^ 

i=l  (l-2t  cosij;  .+t.^^'^^  (5.39) 

i#j  J  J"  J 


because  all  the  terms  are  positive  (til).  Our  goal  is 
to  use  (5.39),  to  determine  a  radius  Rg  of  the  geosphere 
which  will  guarantee  the  convergence  of  the  solution.  It 
would  be  very  difficult  to  use  all  the  terms  of  the  above 
inequality  to  determine  Rg  ,  for  large  N  values.  Below, 
we  present  a  much  simpler 
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method  to  estimate  this  radius,  without  having  to  solve 

a  large  inequality  (such  as  (5.39))  each  time  that  we  seek 
a  solution.  At  first,  note  that  each  term  in  the  sum  above 
is  a  function  of  the  distance  ,  and  of  the  height  of 


(5.9) 


(5.40) 

The  smaller  the  distance  li/..  ,  and  the  smaller  the  parameter 
t  ,  the  larger  the  term  in  tne  sum  becomes.  In  other  words, 
the  maximum  term  in  (5.39)  corresponds  to  'l^miij  ^max  ’ 

Some  typical  values  for  these  terms  are  given  in  table  5.1, 
assuming  R  -  Rb  =  100  meters.  From  these  results  we  can 
see  that  for  very  dense  data  (!'  grid,  or  smaller),  and 
for  high  elevations  (above  4  km),  each  term  in  (5.39)  becomes 
almost  equal  to  1.0,  and  therefore,  the  condition  (5.39) 
might  not  be  satisfied.  However,  (5.39),  as  well  as  its 
original  condition  (5.19),  are  both  sufficient  conditions. 
This  means  that  the  iteration  might  converge,  even  if  this 
condition  is  not  satisfied.  We  will  see  an  example  in  our 
tests  in  the  next  chapter. 

The  Inequality  (5.39)  determines  the  minimum  radius 
of  the  geosphere  (Sg)  which  guarantees  the  convergence  of 
the  iteration.  For  any  radius  Rg  larger  than  this  minimum 
value  Rb  ,  and  such  that  the  geosphere  is  completely  imbed¬ 
ded  inside  the  earth,  the  convergence  is  certain: 

Rb  <  '  guarantees  convergence 

If  all  the  terms  in  (5.39)  had  the  same  magnitude,  we  could 
solve  for  Rb  from  an  inequality  of  the  form: 

(1  -  t.i)^  ^  1 

( 1  -  2tj  cosiir^~  +  t'p5j  “FT 


(5.42) 


tne  point  Oj  k  ana  tnat: 


Rr 

TT+'h' 


which  yields: 


t  . 
min 


ft  +  h 


max 


but  the  magnitude  of  these  terms  varies  very  rapidly,  espec¬ 
ially  with  the  value  of  the  distance  <1;  (see  in  table  5.1). 
As  we  have  already  mentioned  above,  the  largest  value  for 
these  terms  corresponds  to  the  minimum  distance,  and  to 
the  maximum  height. 


Consequently,  if  the  data  is  distributed  on  a  regular 
grid  -  as  in  the  simulations  described  in  the  next  chapter. 
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there  will  be  four  terms  in  (5.39)  whose  magnitudes  are  the 
largest:  these  terms  correspond  to  the  point  Q.  with  the 

maximum  elevation,  and  to  the  four  closests  points  in  the 
vicinity  (figure  5.4). 


Figure  5.4:  The  Point  Qj  with  Maximum  Elevation  h^^^ 
to  which  the  Maximum  Terms  in  (5.39)  Correspond.  ’ 


From  the  tests  reported  in  the  following  chapter,  we 
found  that  it  would  be  sufficient  to  determine  Rg  by  solv¬ 
ing  the  equation  below  (the  term  1/10  was  determined  empir¬ 
ically)  : 


(1  -  t 


mm 


cosil  +  t^r 
^min  min  ^ 


(1  -  2t 


min 


1 

V5 


(5.43) 


where 


'min 


R  +  h 


max 


from  which  it  can  be  easily  found  that: 


t  . 
min 

in  which 


=  P  _ 


1  -  COS;<^n,Hn  (1/10)^^^ 
1  -  (1/10)‘^^ 


(5.44) 


(5.45) 


P  =• 


i  1 


(5.46) 


and  finally: 
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=  (R  + 


^max^ 


[P  -  /P-'-  1  ] 


(5.47) 


This  equation  has  been  used  to  determine  nominal  values 
for  the  depth  (R  -  Rg)  of  the  geosphere,  for  various  values 
of  the  minimum  distance  ’-I'njin  maximum  height  hjiiax' 

R  is  again  the  radius  of  the  mean-earth  sphere,  and  since* 
the  geosphere  must  be  completely  imbedded  inside  the  earth, 
the  iterative  solution  for  the  spikes  will  converge,  if 
we  select  a  radius  Rq  such  that:  Rg  <  Rg  <  R = 6371  km 
(see  condition  (5.41)).  These  nominal  values  for  the  depth 
of  the  minimum  geosphere  are  given  in  table  5.2.  For  those 
blocks  in  this  table  where  no  depth  is  given,  the  data  is 
so  dense,  and  at  such  a  high  elevation,  that  the  iterative 
solution  might  not  converge  (as  a  matter  of  fact,  we  found 
that  in  these  cases  Rg  is  larger  than  R  ) . 

In  the  last  column  of  table  5.2,  we  list  the  correspond¬ 
ing  depth  from  the  Sjoberg's  (1978,  p.64)  "rule",  that  the 
optimum  depth  of  the  geosphere  should  be  half  the  distance 
between  neighboring  surface  points.  We  can  see  that  Sjoberg’s 
depths  agree  very  well  with  ours  (for  the  h  =  0  case).  This 
agreement  will  be  demonstrated  again  in  chapter  9,  where 
some  tests  with  real  data  and  the  Dirac  approach  are  discussed. 
The  advantage  of  equation  (5.47)  in  computing  the  depth 
of  the  geosphere,  is  due  to  the  fact  that  Rb  is  now  a 
function  of  both  the  distance  between  the  observations, 
and  of  the  maximum  elevation  of  the  data  points. 


5.3.3  The  Computation  of  the  Gravity  Disturbance  Vectors 
in  Space  (Dirac  Approach) 

It  remains  to  use  the  spikes  AgS  on  the  geosphere, 
to  compute  the  components  of  at  a  space  point.  This 

is  done  by  substituting  (5.33)  into  (5.24),  and  using  the 
definition  of  the  Delta  Function  (5.34).  Thus  we  obtain 
( note  dS  =  R^  do ) : 


= 

9r  1 

1=1 

-  ^ 

N 

I 

i=l 

.  _s  ,  9S(r,iJ^)  . 

<SSj)p  • 

-  Rb_ 
^*0 

N 

1 

-i 

,_s  ,9S(r,;J;)v 

(5.48) 


30.0  29068  28968  28570  28072  27707  26081  25086  24091  19113  27750 
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where  the  terms  ,  are  given  by  (5.26), 

(5.30),  and  (5. 31) ^respectively .  The  three  components  of 
the  vector  ^  at  P  are  then  computed  from  (5.32).  Finally, 
equations  (3.24)  and  (3.25)  can  be  used  to  compute  tbe  two 
distinct  components  3T/3r  ,  and  (1/r)  dT/d\p  in  the  polar 
(r.  Ip)  coordinate  system. 


5.4  The  Initial  -  Value  Method 


The  mean  value  approach,  and  the  Dirac  approach  which 
have  been  discussed  in  this  chapter  are  the  two  most  rep¬ 
resentative  discretization  techniques  thatcanbe  found  in 
the  literature.  Recently,  a  new  method  has  been  proposed, 
called  the  initial-value  method  (Nakiboglu  and  Lim,  1979). 
From  the  theoretical  standpoint,  this  method  is  similar 
to  the  mean-value  approach,  because  the  given  data  and  the 
reduced  quantities  on  the  geosphere  are  assumed  to  be  given 
as  mean  values  over  a  finite  niunber  of  blocks.  The  number 
(m)  of  the  fictitious  gravity  anomalies  Ag*  on  the  geo¬ 
sphere,  is  not  necessarily  equal  to  the  number  (n)  of  the 
surface  gravity  anomalies  Ag  .  The  solution  for  is 

obtained  from  the  numerical  solution  of  a  set  of  ordinary 
first-order  differential  equations  using  the  iterative  Runge- 
Kutta's  method.  However,  there  is  anumber  of  points  we 
would  like  to  emphasize  concerning  this  approach,  as  compared 
to  the  methods  described  in  this  chapter. 

a.  For  a  real-world  application,  the  Dirac-approach  method 
with  point  anomalies  on  the  surface  of  the  earth,  is 
more  advantageous  because  it  avoids  the  assumptions 

on  the  regular  distribution  of  the  data,  on  the  location 
of  the  points  where  the  kernels  are  evaluated,  and  the 
computation  of  the  areas. 

b.  Nakiboglu  and  Lim  (ibid),  have  restricted  themselves 
in  computations  of  geop-spherop  separations  (C)»  and 
deflections  of  the  vertical  (^  ,  n)  at  the  surface  of 
their  model  only,  and  not  in  space.  Furthermore,  they 

do  say  that  the  Dirac  approach  would  yield  better  results 
than  their  method. 

c.  In  their  numerical  example,  only  a  small  number  of  sub¬ 

divisions  (40  or  80  blocks)  has  been  used  for  the  compu¬ 
tation  of  ?,  5,  and  n  . 

d.  As  we  have  already  shown  in  section  5.3.2,  the  Dirac 
approach  permits  a  relatively  easy  way  for  determining 
the  radius  of  the  geosphere  that  will  guarantee  the 
convergence  of  the  iterative  solution,  Nakiboglu  and 
Lim  (ibid),  kept  the  radius  of  the  geosphere  fixed  at 

1  km  below  the  mean  earth  sphere  R  . 


Chapter  6 

SIMULATION  TESTS  WITH  THE  DIRAC  APPRO\CH 
6.1  Introduction 

In  this  chapter  we  describe  the  simulation  tests  which 
we  performed  in  applying  the  Dirac  approach  to  the  terrain 
model,  just  as  we  did  in  chapter  4  with  the  Green's  approach. 
Similar  simulation  studies  are  reported  in  (Reit,  1966), 
but  they  follow  a  flat-earth  approximation  for  the  5  com¬ 
ponent  of  the  deflection  of  the  vertical  only.  It  should 
be  clear  by  now  that  the  data  needed  for  the  Green's  approach 
is  gravity  disturbances,  and  not  gravity  anomalies  as  in 
the  discrete  Dirac  approach.  Therefore,  the  gravity  distur¬ 
bances  that  are  computed  from  the  model  (equation  3.18), 
have  to  be  transformed  to  anomalies.  This  can  be  achieved 
by  using  the  fundamental  equation  of  physical  geodesy  (see 
in  Heiskanen  and  Moritz,  1967,  p.88): 


the  gravity  disturbance  being  computed  from  (3.18),  and 
the  disturbing  potential  from  (3.16).  This  particular  step 
has  been  ignored  in  similar  simulation  studies  in  the  past 
(Molodensky,  et  al.,  1962,  p.217  ;  Bjerhammar,  1976,  p.43; 
Sjoberg,  1978,  p.37K  all  of  which  use  the  conical  Molodensky's 
model.  However,  the  magnitude  of  the  disturbing  potential 
T  on  the  model  is  not  larger  than  1  .  kgal.m  (see  figures 
3.2  through  3.5),  ana  therefore,  the  magnitude  of  the  second 
term  in  the  right-hand  side  of  (6.1)  is  no  greater  than 
0.3  ragal.  Out  simulation  studies  for  the  Dirac  approach 
consist  of  a  six-step  procedure,  which  can  be  described 
as  follows: 

Step-1 .  The  exact  values  of  the  components  of  ^  are  com¬ 
puted  from  equations  (3.23),  (3.24),  and  (3.25)  at 
selected  points,  over  and  in  the  vicinity  of  the  model. 
Step-2 .  The  gravity  anomalies  are  computed  on  the  surface 
of  the  model  from  equations  (6.1)  and  (3.18).  The 
extent  of  the  area  around  the  model,  and  the  grid  spacing 
are  specified,  and  the  data  is  generated  on  a  regular 
grid. 

Step-3 .  The  radius  of  the  (smallest)  geosphere  that  will 
guarantee  the  convergence  of  the  iterative  solution 
is  computed  from  (5.47). 
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Step-4 .  The  Gauss-Seidel  iterative  method  is  applied  for 
the  computation  of  the  gravity  spikes  AgS  at  the 
carrier  points. 

Step-5.  The  components  (6gj.  ,  6g_  ,  <5gx)  are  obtained 

from  (5.48),  and  them,  (5, '»2)^yields  the  three  compon¬ 
ents  8T/3X  ,  3T/3Y  ,  3T/3Z  at  the  selected  space 

points.  Then,  (3.24)  and  (3.25)  can  be  used  to  compute 
the  two  distinct  components  3T)3r  ,  (1/r)  3T/3\iJ  in 

the  polar  (r  ,ii/)  system. 

Step-5*  (optional).  The  N  spikes  on  the  geosphere  must 
satisfy  all  the  surface  gravity  anomalies.  Equation 
(5.35)  can  be  used  to  compute  the  surface  anomalies 
from  the  spikes,  which  can  be  compared  ot  the  original 
data.  Any  difference  between  the  given  data  and  the 
computed  anomalies  are  due  to  the  fact  that  the  iteration 
has  been  forced  to  stop  before  a  16-digit  agreement 
was  reached  between  successive  approximations. 

Step-6 .  The  components  from  steps  1  and  5  are  compared 

in  terms  of  their  relative  percentage  difference  (from 
equations  4.1)). 


6.2  Simulation  Tests 


The  difficulty  in  applying  the  Dirac  approach  arises 
from  the  fact  that  an  iterative  procedure  is  used  to  solve 
for  the  gravity  spikes  on  the  geosphere  (section  5.3.1). 

In  terms  of  computer  time  (on  an  Amdahl  470  machine),  it 
takes  about  16  seconds  to  complete  just  one  iteration  for 
576  spikes.  Therefore,  it  would  be  practically  impossible 
to  handle  a  large  amount  of  surface  data  (say  more  than 
2000  data  points),  because  of  the  large  CPU  time  that  is 
required  to  solve  such  a  large  system  of  equations.  This 
is  why  we  decided  to  make  the  simulation  tests  for  the  Dirac 
approach,  using  a  relatively  small  amount  of  data  points 
(up  to  576  surface  gravity  anomalies).  The  characteristics 
of  the  models  are  defined  from  the  parameters  listed  in 
table  3.1.  A  description  of  the  tests  is  given  in  table 
6.1  for  easy  comparison.  It  must  be  kept  in  mind  that  the 
iterative  solution  for  the  gravity  spikes  can  be  carried- 
out  to  any  desirable  level  of  accuracy,  which  is  controlled 
by  the  convergence  criterion  (5.38),  and  by  the  maximum 
number  of  iterations  allowed  (kj^g^j^). 

The  question  is  now  how  many  iterations  are  actually 
needed  to  compute  the  gravity  spikes.  In  other  words,  a 
12-digit  agreement  between  successive  iterations  at  all 
spikes  according  to  conditon  (5.38)  might  be  too  much,  and 
that  the  iteration  should  rather  be  terminated  whenever 
the  RMS  difference  between  the  surface  anomalies,  and  the 
anomalies  computed  from  the  spikes  at  the  same  points,  is 
below  a  certain  level  (say  1.  mgal ) . 
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In  order  to  investigate  the  effect  of  the  number  of  iter¬ 
ations,  on  the  accuracy  of  the  components  of  o  ,  three  tests 
were  made  using  the  spherical  cap  model  (576  anomalies  in 
a  0.4°x 0.4°  area) ,  and  terminating  the  iteration  process 
after  8,  10,  and  12  iterations  respectively.  The  exact  com¬ 
ponents  from  the  model,  and  the  components  computed  from 
the  spikes  at  7  space  points  (at  5  km  altitude  above  the 
mean  sphere  R  ),  are  listed  in  table  6.2.  The  differences 
between  the  exact  and  the  computed  component  ,  are  also  listed 
in  parenthesis.  The  last  two  columns  of  table  6.2  contain: 

(a) .  The  RMS  difference  between  the  exact  components,  and 

the  computed  ones,  and 

(b) .  The  RMS  difference  between  the  original  surface  anomalies 

and  those  computed  at  the  same  points  from  the  spikes, 
after  the  iteration  is  terminated. 


The  fact  that  these  RMS  differences  become  smaller  as 
the  number  of  iterations  increases,  indicates  that  the  num¬ 
erical  solution  for  the  gravity  spikes  probably  converges 
(the  relative  percentage  errors  from  the  kn,ax“^2,are  given 
in  table  6.4,  and  as  we  can  see  are  quite  small).  Therefore, 
we  decided  to  terminate  the  iteration  for  all  tests  at  kniax= 
12,  for  which  the  RMS  differences  between  the  original  and 
the  computed  anomalies  are  smaller  than  1.  mgal  (see  also 
table  6.1). 

The  radius  of  the  geosphere  is  computed  from  equation 
(5.47),  but  for  most  of  our  tests  with  very  dense  data  (1' 
grid,  see  table  6.1),  R5  was  forced  to  be  50  meters  smaller 
than  R  .  This  had  to  be  done,  because  for  very  dense  data 
at  high  elevations  (as  on  our  4.1-km  models),  the  radius 
Rg  from  (5.47)  is  larger  than  R  .  (note:  the  50-meters 
depth  for  the  geosphere  was  selected  on  a  rather  arbitrary 
basis,  since  the  geosphere  must  be  completely  imbedded  inside 
the  earth,  and  therefore  inside  the  mean-earth  sphere  R  ). 

This  paradoxial  result,  i.e.  Rg  larger  than  R  ,  actually 
means  that  this  kind  of  data  does  not  satisfy  the  sufficient 
condition  (5.39),  and  hence,  the  iterative  solution  for  the 
spikes  might  not  converge.  However,  the  tests  that  follow, 
’ndicate  that  even  in  this  case  the  iteration  converges, 
which  is  explained  by  the  fact  that  (5.39)  is  only  a  sufficient 
condition  for  convergence. 

Let  us  first  start  the  discussion  on  the  simulations 
with  the  errors  from  the  spherical  cap  model,  already  used 
above  for  the  determination  of  the  maximum  number  of  iter¬ 
ations.  A  total  number  of  576  anomalies  were  generated  on 
its  surface,  within  a  2°  x  2°  area  on  a  5'  grid.  This  data 
implies  a  radius  Rg  753  meters  smaller  than  R  .  The  errors 
in  the  components  of  0  at  5  km  altitude  are  large  (table 
6.3),  but  the  reason  is  obvious:  most  of  the  gravity 
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information  generated  by  the  model  is  missing,  since  the 
data  points  are  too  far  away  from  it.  As  we  can  see  in  table 
6.4,  by  using  data  on  a  1*  grid  within  a  0.4°  x  0.4°  area  - 
576  points  again-,  the  errors  are  twc  orders  of  magnitude 
smaller  than  with  data  on  a  5*  grid  in  a  much  larger  area. 

The  errors  from  these  two  tables  indicate  that  for  the  very 
smooth  topography  of  the  spherical  cap  model,  the  results 
from  the  Dirac  approach  are  very  sensitive  to  the  density 
of  the  surface  data,  but  the  method  works  very  well  with 
dense  data  (say  on  a  1 '  grid),  provided  that  this  data  is 
not  too  dense  if  we  want  to  avoid  divergence  of  the  iter- 
at ion . 

For  a  10°-inclination  conical  model,  with  576  anomalies 
ona  1'  grid  as  above ,  the  erro^s  at  5,  10,  20,  and  100  km 
altitude  are  given  in  tables  6.5,  6.6,  6.7,  and  6.8  respec¬ 
tively.  At  altitude  5  ^m,  the  errors  from  the  spherical 
cap  model  are  slightly  smaller  than  the  errors  from  the  10°- 
inclination  conical  model  (table  6.4  vs.  6.5).  It  is  quite 
remarkable  that  even  over  the  edge  of  the  area  at  5  km  alti¬ 
tude,  the  errors  are  no  more  than  9%.  Not  only  that,  but 
comparing  with  the  results  from  the  Green's  approach  at  10 
km  are  almost  three  times  smaller,  despite  the  fact  that 
the  amount  of  data  is  now  two  orders  of  magnitude  smal¬ 

ler.  More  remarkable  is  the  fact  that  directly  above  the 
model,  at  5  km,  10  km,  and  20  km  altitude,  the  errors  are 
less  than  3",  as  opposed  to  almost  25^  errors  from  the  Green's 
approach.  However,  as  the  altitude  increases,  the  errors 
from  the  Dirac  approach  become  slightly  larger,  in  contrast 
to  the  results  from  the  Green's  approach.  Nevertheless, 
at  100  km  the  errors  from  the  Dirac  approach  are  almost  four 
times  smaller  than  the  errors  from  the  Green's  approach  (table 
6.8  vs .  4.16). 

An  overall  examination  of  \,he  two  techniques,  using 
identical  data  and  models  (10°-cone,  0.4°  x  0.4°  area.  1' 
grid,  576  anomalies),  shows  that  the  Dirac  approach  in  general 
is  superior  to  the  Green's  approach  in  terms  of  the  magnitude 
of  the  errors.  The  problem  with  the  Dirac  approach  is  that 
it  will  fail  with  very  dense  data  at  high  elevations,  ard 
it  requires  considerably  more  computer  time  because  of  the 
iterative  solution  for  the  gravity  spikes. 

For  a  40° - inc 1 inat ion  conical  model  (table  6.9),  the 
errors  are  quite  large,  reaching  almost  90%  over  the  whole 
area.  However,  we  have  to  realize  that  the  extent  of  a  40° 
cone  with  4.1  km  height  is  =  2.64',  and  thus,  there  are 
only  16  data  points  on  the  conical  surface  with  a  1'  grid. 
Apparently,  this  amoung  of  data  is  not  enough  to  represent 
the  model's  gravity  field  (figure  3.4)  for  the  Dirac  approach. 
On  the  other  hand,  the  Green's  approach  with  the  same  model 
and  data,  and  at  the  same  altitude,  has  yielded  considerably 
smaller  errors  (table  4.11).  Attempts  to  use  smaller  grid 
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spacing  with  the  Dirac  approach  failed,  because  of  the  diver¬ 
gence  of  the  iterative  solution  for  the  spikes. 

All  the  tests  in  the  present  study  use  synthetic  data 
on  the  surface  of  the  models,  and  therefore,  there  are  no 
data  errors  (the  gravity  anomalies  are  rigorously  computed 
from  the  two  point  masses).  In  addition,  and  similarily 
to  the  simulations  with  the  Green's  approach,  due  to  the 
very  local  characteristics  of  the  models( f ip  are  3.2  through 
3.5),  there  are  no  significant  truncation  errors  (i.e.  errors 
caused  by  neglecting  the  information  outside  of  the  working 
area).  The  resulting  errors  which  are  given  in  tables  6.3 
through  6.9  are  due  only  to  the  following  factors: 

(a) .  Errors  in  the  numerical  solution  for  the  gravity  spikes, 

and 

(b) .  Errors  caused  by  using  a  small  number  of  data  points 

on  the  surface  of  the  model. 


From  our  experience  with  the  simulations  above,  we  can 

draw  the  following  conclusions: 

1.  As  the  altitude  of  the  space  point  increases,  the  errors 
increase  too.  This  effect  is  opposite  to  that  from  the 
Green's  approach.  Possibly,  the  small  amount  of  data 
points  on  the  surface  of  the  models  (determined  by  the 
grid  spacing),  is  responsible  for  this  discrepancy. 

2.  In  general,  the  errors  over  the  model  seem  to  be  smaller 
than  at  the  distant  points,  again  in  contrast  to  the 
Green's  approach. 

3.  The  finer  the  grid  spacing,  the  smaller  the  errors  become, 
but  the  disadvantage  of  the  Dirac  approach  is  that  it 
will  not  converge  with  very  dense  data  at  high  elevations. 

4.  The  larger  the  inclination  of  the  model,  the  larger  the 
errors  become,  again  in  contrast  to  the  results  from 
the  Green's  approach.  It  seems  that  the  Dirac  approach 
is  more  sensitive  to  high  inclinations  than  the  Green's 
approach.  For  a  smooth  topography  (the  case  of  the  spher¬ 
ical  cap  model,  and  of  the  10®-inclination  cone),  the 
Dirac  approach  is  superior  to  the  Green's  approach. 
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Table  6.5;  Comparison  of  Gravity  Disturbance  Vector  Components  in  Space 
Computed  Rigorously,  and  From  the  Dirac  Discrete  Approach 
(values  in  mgal ;  height  in  meters) 
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Table  6.7:  Comparison  of  Gravity  Disturbance  Vector  Components  in  Space 
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6.3  Comparison  of  the  Green’s  with  the  Dirac  Approach 

In  section  4.2  and  6.2  we  prsented  the  simulation  studies 
with  the  two  approaches.  Let  us  now  review  the  results, 
and  compared  the  two  methods. 

The  Green's  approach  requires  as  data  gravity  distur¬ 
bances  dg  on  the  surface  of  the  earth,  the  disturbing  poten¬ 
tial  T  ,  the  elevation  h  ,  and  the  two  components  of  the 
inclination  at  these  points.  Actually,  SS  >  T  ,  and  h 
are  referred  to  the  center  of  the  blocks,  after  a  subdivision 
of  the  surface  has  been  made  using  a  certain  grid  spacing, 
and  the  whole  approach  is  thus  an  approximation.  The  fact 
that  the  errors  of  J  are  not  reduced  by  using  a  finer  grid, 
is  an  indication  that  this  kind  of  approximation  (combined 
with  the  neglect  of  2nd  order  variations  of  the  topography) 
is  not  sufficient .  The  method  does  not  work  on  the  surface, 
and  close  to  it  gives  very  large  errors.  As  the  altitude 
of  the  space  point  increases,  the  errors  descrease.  Because 
of  the  very  local  characteristic  of  the  model  used  in  the 
simulations,  the  increase  of  the  integration  area  doesnot 
reduce  the  errors,  and  hence,  there  are  no  significant  trunca 
tion  errors.  Also,  as  the  inclination  of  the  model  increases 
(beyond  20°),  the  errors  become  smaller.  The  method  is  much 
faster  than  the  Dirac  approach,  but  it  requires  the  data 
be  given  on  a  regular  grid. 

The  Dirac  approach  requires  as  data  only  surface  gravity 
anomalies  and  elevations  at  points,  not  necessarily  on  a 
regular  grid.  The  method  requires  an  iterative  solution 
for  the  anomalies  on  the  geosphere,  which  might  diverge  for 
dense  data  at  high  elevations.  As  the  inclination  of  the 
model  increases,  the  errors  increase.  The  same  happens  when 
the  altitude  of  the  space  point  increases,  but  the  errors 
are  stil  smaller  than  those  from  the  Green's  approach.  In 
general,  this  method  works  well,  both  close  to  the  surface 
of  the  earth,  and  at  high  altitudes.  Its  basic  disadvantage 
is  the  fact  that  it  cannot  be  used  for  large  amounts  of  data, 
since  it  involves  an  analytical  solution  for  the  anomalies 
on  the  geosphere,  which  is  a  very  time-consuming  process. 
Also,  it  might  diverge  for  very  dense  data  at  high  elevations 

If  it  is  a  matter  of  choice  between  the  two  methods, 
from  our  experience  with  these  simulation  studies,  we  would 
recommend  the  Green's  approach  for  high  altitudes,  with  a 
large  number  of  data  on  the  surface  of  the  earth,  and  the 
Dirac  approach  close  to  the  surface,  if  a  small  number  of 
point  data  is  available. 


6.4  Related  Work 
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A  number  of  authors  have  applied  Bjerhammar's  methods 
for  the  computation  of  the  external  gravity  field  of  the 
earth.  Some  of  them  followed  the  original  method  (based 
on  the  mean-value  approach),  and  some  other  have  developed 
methods  based  on  the  fitting  of  the  surface  dr.ta  to  an  m- 
degree  polynomial.  To  the  best  of  my  knowledge,  the  Dirac 
approach  has  never  been  applied  to  real  or  simulated  data 
for  the  computaion  of  the  complete  gravity  vector  in  space, 
without  neglecting  the  topography. 

The  earliest  work  on  this  subject  was  published  by  Reit 
(1966).  It  is  a  flat-earth  (planar  approximation)  simulation 
study,  using  two  models,  one  of  which  is  almost  identical 
to  the  conical  model  described  in  chapter  3.  The  gravity 
anomalies  on  the  geosphere  ( Ag* ) ,  are  computed  iteratively, 
and  the  C  -component  of  the  deflection  of  the  vertical 
is  computed  and  compared  to  its  exact  value  on  the  model's 
surface,  and  in  space.  The  errors  are  almost  two  times  larger 
than  those  found  from  our  analysis.  For  example,  from  Reit's 
tables  5  and  6,  the  errors  at  5  km  above  the  cone  are  almost 
10%,  while  at  the  same  altitude  the  errors  from  the  Dirac 
approach  are  less  than  2%  (tables  6.4  and  6.5).  For  reasons 
that  will  be  explained  in  the  chapter  9  we  have  to  point 
out  that  Reit  computed  the  C  -component  on  the  surface  of 
the  model  at  the  same  points  where  the  "synthetic"  data  is 
given. 

Forstner  (1966)  used  125  gravity  anomalies  in  Cyprus, 
and  performed  four  polynomial  fittings  to  this  data.  However, 
these  methods  were  applied  in  computing  surface  components 
of  the  gravity  vector  only.  In  addition,  the  RMS  difference 
between  the  original  data,  and  the  anomalies  computed  from 
the  reduced  (Ag*)  data,  was  found  to  be  between  7  and  22 
mgal  (Forstner,  ibid,  p.55).  This  difference  is  much  larger 
than  that  found  from  our  tests  (see  table  6.1),  despite  of 
the  fact  that  the  maximum  elevation  on  our  model  is  4.1  km, 
as  opposed  to  1.8  km  over  Cyprus. 

Barlik  (1971)  describes  other  methods  which  are  proposed 
for  areas  with  high  elevations  in  order  to  avoid  the  divergence 
problems  which  are  encountered  by  Forstner  (ibid).  An  auxil- 
liary  sphere  is  used,  not  imbedded  inside  the  earth  as  the 
geosphere,  that  the  terms  t ^  = Rg/r j  are  not  very  small 
for  the  points  at  high  elevations.  However,  the  proposed 
method  is  not  fully  automatted,  and  the  evaluation  of  terrain 
corrections  is  required. 

The  original  work  by  BJerhammar  (1963,  1975,  1976,  1978), 
and  that  by  Sjoberg  (1975,  1978),  have  already  been  mentioned 
in  section  5.1.  BJerhammar  (1976)  has  presented  some  results 
from  the  application  of  the  Dirac  approach  for  the  computation 


of  deflections  of  the  vertical  on  the  surface  of  the  earth. 
Finally,  Sjoberg  <1975,  1978)  has  worked  with  a  limited  amount 
of  real  data  for  the  computation  of  gravity  anomalies  on 
the  surface  of  the  earth.  Because  we  feel  that  these  results 
need  some  additional  elaboration,  we  will  repeat  his  tests 
in  chapter  9  with  some  comments. 


Chapter  7 


COMPUTATION  OF  THE  GRAVITY  DISTURBANCE  VECTOR  IN 
SPACE  USING  THE  METHOD  OF  LEAST  SQAUARES  COLLOCATION 

7.1  Introduction 


Let  us  now  apply  the  method  of  I east-Squares  Collocation 
to  the  estimation  of  the  components  of  ^  ,  from  the  simu¬ 
lated  data  on  the  surface  of  the  nrodel.  As  with  the  other 
two  methods  in  this  report,  it  is  easy  to  evaluate  the  appli¬ 
cability  of  collocation  to  our  problem,  because  the  exact 
components  of  t  can  be  computed  directly  from  the  model. 

A  great  number  of  papers  have  been  published  on  appli¬ 
cations  of  collocation  to  the  determination  of  the  gravity 
field  of  the  earth,  on  its  surface  as  well  as  in  the  exterior 
space.  For  the  estimation  of  gravity  anomalies,  geoidal 
undulations,  and  deflections  of  the  vertical  see  for  example 
in  (Lachapelle,  1977;  Sjdberg,  1978;  Tscherning  and  Fors- 
berg,  1978;  Forsberg  and  Tscherning,  1981).  For  computation 
of  gravity  anomalies  in  space  see  Rapp  and  Hajela  (1975). 

For  the  computation  of  the  covariance  functions  -  the  most 
essential  quantities  in  any  estimation  through  collocation, 
see  Tscherning  and  Rapp  (1974),  Tscherning  (1976),  and  Siinkel 
(1979).  Advanced  aspects  on  collocation  are  discussed  in 
Moritz  (1980),  and  Moritz  and  Siinkel  (1978). 

The  collocation  estimation  of  a  quantity  s  from  gravity 
anomalies  Ag  ,  can  be  done  using  the  equation  below  (Moritz, 
1972) : 


s 


where : 


C  *  (C.  .  +  D)-^  Ag 

-s,Ag  -Ag,Ag 


(7.1) 


“Agf Ag 
D 


-s,Ag 


the  vector  of  the  m  observations  of  gravity 
anomalies 

the  covariance  matrix  (mxm)  of  the  observed 
anomalies. 

the  covariance  matrix  (mxm)  of  the  measuring 
errors  in  the  anomalies.  This  is  taken  to  be  a 
null  matrix,  since  we  assume  errorless  observa¬ 
tions  for  our  simulation  study, 
the  cross-covariance  vector  (m),  between  the 
quantity  s  to  be  predicted,  and  the  gravity 
observatibns. 
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The  predicted  s  might  be  any  quantity  as  long  as  the 
corresponding  covariance  function  between  s  and  the  data 
(Ag)  exists.  For  example,  s  may  be  geoidal  undulation, 
gravity  anomaly,  first  or  higher  order  derivative  of  the 
disturbing  potential,  etc.  The  two  problems  that  we  encounter 
in  any  collocation  prediction  are  (a)  the  inversion  of 
a  large  matrix  IC  +  D)  in  case  of  a  large  number  of  data, 
and  (b)  the  computation  of  the  auto  -  and  cross-covariances 
between  the  various  quantities  in  (7.1). 

For  the  computation  of  the  various  covariances  we  decided 
to  use  the  covariance  approximation  procedure  which  is  doc¬ 
umented  in  Sunkel  (1979).  In  fact,  we  have  used  a  new  version 
of  his  computer  program,  documented  in  Sunkel  (1980),  the 
difference  between  the  old  and  the  new  version  being  minor. 

The  whole  procedure  is  based  on  the  generation  of  a  suffic¬ 
iently  dense  network  of  covariances,  and  then,  on  the  computa¬ 
tion  of  any  other  covariance  at  a  point  inside  this  network 
by  a  differentiation  -  interpolation  procedure  with  spline 
functions.  The  accuracy  of  this  procedure  can  be  as  high 
as  we  wish,  depending  on  the  density  or  the  covariance 
network.  Interpolating  at  any  other  point  from  the  grid 
values,  rather  than  computing  the  covariance  directly,  re¬ 
duces  the  CPU  time  significantly  (Sunkel,  1979,  p.22). 


7.2  Collocation  Prediction  of  the  Gravity  Disturbance  Vector 

The  computer  program  mentioned  above,  computes  the 
auto-  and  cross-covariances  between  fourteen  quantities, 
among  which  we  find  the  radial  component  of  the  disturbing 
potential,  and  the  components  of  the  deflection  of  the  verti¬ 
cal  : 

3T 

“  — in  Ebtvbs  (1  E  =  10”®sec”^) 

i  in  arcsec 

n  in  arcsec 

Ag  in  mgal 

For  simplicity,  let  us  denote  by  g  the  product  C~^  ^ 
in  (7.1); 

cov(Agi,Agi)  .  .  .  cov(Agi,Ag  )  Ag 

#  ro  I 

I  f 

cov(Agjjj,Agi)  .  .  .  cov(Agjjj,Agjjj) 

(7.2) 


95. 

Then,  the  application  of  (7.1)  to  the  determination 
of  the  quantities:  -OT/3r)/r  ,  ^  ,  and  n  yields: 

-(l/r)aT/3r  =  [cov(-l/r)3T/3r,  Agi],...]  Q  ,  in  Eotvos 
C  =  [cov(C,Agi),  .  .  .}  Q  ,  in  arcsec 


n  -  [cov(n,Agi),  .  .  .]  g  ,  in  arcsec 

The  partial  derivatives  of  the  disturbing  potential  with 
respect  to  r  ,  cf  ,  and  X  ,  are  then  given  from  the  well- 
known  equations  (Siinkel,  ibid,  p.21): 


3T 

3<{) 

3T 

3X 


-r(-LLilL) 

r 

-r  Y  5 

-ry  cosi  n 


,  in  10~®m.sec“^ 

,  in  m  .sec“^  .  arsec 


,  in  m  .sec~^.  arsec 

(7.4) 


The  components  of  T  can  now  be  transformed  to  compon¬ 
ents  of  the  gravity  disturbance  vector  ^  ,  using  (5.22): 


6g  .  |JL  =  -r  10- 

v&r  9  j.  ^  r 


,  in  m.sec 


,  in  m.sec"^ 


,  in  m.sec  ^ 

(7.5) 


Then,  equation  (5.32)  is  used  to  rotate  the  components 
of  to  a  geocentric  (X  ,  Y  ,  Z)  system: 


3T 

5gx 

cosAp 

cosAp 

3X 

cos<i>p 

-sin^p 

-sinAp 

«gr 

3T 

<5gY 

sin  Ap 

~5T 

COS^p 

-sin5 

P 

sinAp 

cosAp 

3T 

~1T 

6gz 

sin^p 

cosip 

0 

7.6l 
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where  the  index  P  means  that  the  computations  are  made 
at  the  space  point  P  (rp  ,  ifp  ,  Xp). 

riuall^i  aiitl  tail  Lc  i4oct>l  i>> 

rotate  these  components  to  the  polar  (r  ,  ij/)  coordinate 
system  in  order  to  obtain  the  two  distinct  components  9T/3r  , 
and  (1/r)  for  the  simulation  tests. 


7.3  Simulation  Tests  with  the  Collocation  Approach 

As  we  described  in  sections  4.2  and  6.2  for  the  Green's 
and  for  the  Dirac  approaches  respectively,  the  terrain  model 
of  chapter  3  is  used  again,  and  the  gravity  anomalies  are 
generated  on  its  surface  from  the  two  point  masses  to  be, 
used  in  the  collocation  prediction. 

The  covariance  network  which  is  the  most  time-consuming 
part  in  the  covariance  approximation  procedure  just  mentioned, 
is  also  generated  in  space,  from  the  surface  of  the  sphere 
up  to  an  altitude  of  111  km  in  radial  direction,  and  from 
the  center  of  the  model  up  to  125’  in  spherical  distance. 

The  grid  has  to  be  denser  for  smaller  altitudes  than  for 
higher  altitudes,  and  denser  for  smaller  spherical  distances 
than  for  larger  ones  (S'unkel,  1979,  p.5).  From  our  prelimin¬ 
ary  tests  with  the  collocation  approach,  we  found  that 
a  non-uniform  grid  as  it  is  shown  in  figure  7.1  would  be 
sufficient  for  the  prediction  of  the  components  of  the  grav¬ 
ity  disturbance  vector  up  to  an  altitude  of  100  km.  The 
grid  that  we  constructed  has  a  spacing  of  1  km  up  to  the 
11  km  altitude,  and  then  the  spacing  becomes  10  km  up  to 
the  maximum  111  km  radial  distance.  In  the  spherical  direc¬ 
tion,  the  grid  spacing  is  O'. 5  for  the  first  5',  then  it 
is  1'  for  the  next  20',  and  finally  it  becomes  10'  up  uu 
the  maximum  spherical  distance  of  125'.  It  takes  approxi¬ 
mately  18  CPU-seconds  in  an  Amdahl  470/6-II  maching  for 
the  construcion  of  such  a  grid. 

The  Tscherning  and  Rapp  (1974  ,  4th  model)  anomaly 
degree  variance  model  was  used  for  the  computation  of  the 
covariance  network.  This  model  corresponds  to  the  2nd  model 
in  Sunkel's  program  (see  also  equation  (9.7),  section  9.3). 
Preliminary  tests  indicated  that  the  use  of  a  "global"  versus 
a  "local  "  nt'i-degree  covariance  function  (the  first  n  anom¬ 
aly  degree  variances  are  set  equal  to  zero),  is  not  a  critical 
factor  for  the  computation  of  t  •  The  use  of  a  global, 
of  a  20th-degree,  and  of  a  40th-degree  covariance  functions, 
resulted  to  components  of  J  ,  which  are  different  by  no 
more  than  0.5  mgal.  The  fact  that  the  collocation  prediction 
is  not  sensitive  to  the  covariance  function  being  used, 
has  been  verified  in  numerous  applications  of  collocation 
(see  for  example  in:  Rapp  andAgajelu,  1975,  p.9;  Lachapelle, 
1977;  Rapp,  1979-a;  Katsambalos,  1980). 
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One  might  argue  that  the  covariance  function  being  used 
should  be  consistent  with  the  synthetic  data  on  the  models. 
Theoretically,  this  is  necessary  due  to  the  fact  that  the 
center  of  mass  of  the  models  does  not  coincide  with  the 
center  of  the  mass  of  the  earth  (the  center  of  the  sphere 
R  in  our  case).  The  construction  of  a  covariance  function 
is  a  very  laborious  procedure,  especially  for  our  simulation 
tests  where  a  great  number  of  different  models  are  used. 

In  addition,  as  we  mentioned  above,  from  ou .  preliminary 
tests  it  was  found  that  the  predicted  quantities  are  not 
sensitive  to  the  covariance  function  being  used.  Therefore, 
for  simplicity  we  decided  to  use  the  global  Tscherning- 
Rapp  covariance  function. 

After  the  covariance  network  is  set-up,  the  equations 
of  the  previous  section  can  be  used  for  the  collocation 
prediction  of  the  components  of  ^  .  The  radii  rp  ,  Tq 
are  computed  as  rp  =  R  +  hp  ,  r^  =  R  +  h^  ,  where  hp  ,  hg 
are  the  heights  of  the  space  pornt  P  ^of  the  surface  point 
Q  above  the  mean-earth  sphere  R  . 

In  table  7.1  the  relative  percentage  errors  of  the 
vector  ^  at  5  km  altitude  are  shown  in  the  usual  form. 

256  anomalies  are  generated  on  the  surface  of  a  10°  conical 
model,  on  a  3'  grid  within  a  0.8°x  0. 8°  area.  These  errors 
are  about  20%  above  the  cone,  they  reach  a  maximum  of  70%, 
and  then  they  fall  below  the  10%  level.  By  decreasing  the 
grid  spacing  to  2'  within  the  same  area  (576  anomalies  on 
the  surface  of  the  model  now),  the  errors  at  5  km  are  now 
10%  (table  7.2),  but  they  grow  larger  as  we  go  away  from 
the  model's  axis  towards  the  boundary  of  the  area. 

Table  7.3  shows  the  errors  at  10  km  altitude  using 
the  10°  conical  model  again.  Comparing  these  errors  with 
those  at  5  km  (table  7.1),  we  see  that  they  are  almost  an 
order  of  magnitude  smaller  right  above  the  model’s  center, 
but  they  are  larger  away  from  it.  At  higher  altitudes  (100 
km  and  above)  the  errors  from  the  collocation  approach  were 
found  to  be  very  large  (more  than  100%).  Any  attempt  to 
increase  the  size  of  the  area  with  2'  or  3'  grid  spacing, 
would  require  the  inversion  of  large  matrix  C  .  Not  that 
it  takes  about  4  minutes  CPU  time  to  form  and  invert  a  576  X 
576  symmetric  matrix  of  the  covariances  (the  IMSL  subroutine 
LINV3P  was  used). 

For  a  direct  comparison  between  the  Dirac  and  the  col¬ 
location  approaches,  the  components  of  t  at  5  km  altitude 
are  tabulated  in  table  7.7  using  the  same  model  (10°  cone), 
and  the  same  data  (576  anomalies  on  a  2’  grid,  within  a 
0 . 8°x  0 . 8°  area) .  The  results  from  the  collocation  approach 
are  in  a  better  agreement  with  the  exact  components  of  t 
right  above  the  model's  center,  but  the  Dirac  approach  is 
obviously  superior  at  the  other  points. 


ns. 

The  errors  from  a  40®-inclination  model  (table  7.4  , 
compared  to  those  from  a  10®-model,  (table  7.2),  are  larger 
over  the  cone,  but  smaller  away  from  it.  In  order  to  dec reap 
the  errors  at  the  points  away  from  the  model's  axis,  we 
should  increase  the  working  area  beyond  l“xl“,but  that 
would  require  the  inversion  of  a  very  large  matrix.  There¬ 
fore,  the  application  of  collocation  to  our  problem  seems 
to  be  a  very  time  consuming  approach. 

For  the  very  smooth  topography  of  the  spherical  cap 
model  which  was  described  in  section  3.^  the  errors  are 
smaller  when  compared  to  those  from  the  conical  model,  and 
they  become  almost  zero  at  point  away  from  the  model  (tables 
7.5  and  7.6).  This  is  true  not  only  at  5  km  altitude,  but 
also  at  10  km  altitude. 

In  order  to  investigate  the  effect  of  the  topography 
on  the  collocation  approach,  we  made  an  additional  test 
for  which  the  elevations  of  the  surface  points  were  forced 
to  be  equal  to  zero(the  remaining  specifications  for  this 
test  are  identical  to  those  in  table  7.2).  The  results 
from  this  test  are  given  in  table  7.8.  Comparing  the  errors 
given  in  tables  7.2  and  7.8,  we  see  that  the  effect  of  the 
topography  is  very  significant,  since  the  errors  from  the 
new  test  are  much  larger  above  the  model's  center,  than 
those  in  table  7.2. 

From  our  experience  from  the  simulation  tests  with 
the  collocation  prediction  of  ,  we  can  make  the  following 
conclusions. 

(a) .  Collocation  is  a  very  time  consuming  approach,  espec¬ 

ially  in  handling  a  large  amount  of  surface  data. 

(b) .  The  errors  in  ^  are  very  large  for  altitudes  larger 

than  10  km  (relative  percentage  errors  at  the  100% 
level  and  above) 

(c) .  For  data  given  on  a  very  smooth  topography,  the  errors 

at  5  km,  and  10  km  altitude  are  relatively  small  (com¬ 
pared  to  the  errors  from  the  conical  models),  but 
not  as  small  as  those  from  the  Dirac  approach. 

(d) .  On  an  overall  basis,  the  Dirac  approach  seems  to  be 

superior  as  compared  to  the  collocation  approach. 
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Table  7.3;  Comparison  of  Gravity  Disturbance  Vector  Components  in  Space 
Computed  uigorously,  and  From  the  Collocation  Approach 
(values  in  mgal;  height  in  meters) 
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Chapter  8 


COMPARISON  OF  THE  IMPROVED  TECHNIQUES 
WITH  THE  CLASSICAL  APPROACH 

As  we  have  already  mentioned  in  chapter  1  ,  there  exist 
three  methods  for  the  computation  of  t  (using  gravity 
data  reduced  to  the  geoid),  which  were  called  the  "classical" 
approaches : 

1.  the  Direct  Integration  Method, 

2.  the  Coating  Me'-'iod,  and 

3.  the  Upward  Cot  t:.’uation  Method. 

The  theory  behind  these  methods  can  be  found  in  Hirvonen 
and  Moritz  (1963"+;  Heiskanen  and  Moritz  (1967);  Mueller 
(1966).  The  Direct  Integration  method  is  computationally 
the  most  difficult  to  use,  but  it  requires  the  least  amount 
of  data,  since  the  other  two  methods  require  -  in  addition 
to  the  gravity  anomalies  which  are  needed  by  all  three  -  , 
geoidal  undulations,  and  deflections  of  the  vertical  on 
the  physical  surface  of  the  earth.  A  computer  program  which 
is  based  on  the  Direct  Integration  Method  is  documented 
in  (Rapp,  1966). 

We  decided  to  use  the  Direct  Integration  Method  as 
our  "classical"  approach,  because  of  its  simplicity  in  terms 
of  data  required.  The  corresponding  equations  for  the  three 
components  of  J  ,  at  a  point  P  (^p  X  p  ,  rp)  in  space 
are  (Heiskanen  and  Moritz,  1967,  p.234): 
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(8.1) 

Strictly  speaking,  the  gravity  anomalies  Ag  =  g^  -  vg 
in  the  equations  above  refer  to  the  geoid;  gQ  is  the  gravity 
reduced  in  free-air  from  the  surface  of  the  earth  to  the 
geoid,  and  yk  the  normal  gravity  on  the  ellipsoid. 

It  must  be  pointed-out  that  in  applications  with  real  data, 
the  normal  gradient  dy/dr  »  -0.3086  ragal/meter  is  used 
for  this  type  of  reduction,  instead  of  the  gradient  3g/3r 
Therefore,  the  resulting  anomalies  are  approximately  surface 
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free-air  gravity  anomalies.  However,  the  difference  between 
the  surface  anomaly,  and  the  anomaly  at  the  geoid  is  small 
Mhifl,  p.Piil).  ronsr'qiif'nt  1 V .  for  our  simulation  tests  with 
the  classical  approach  we  will  be  using  gravity  anomalies 
generated  by  the  model  on  its  surface  (see  equation  (6.1)) , 
similarily  to  the  computations  performed  by  Molodensky  et 
al.  (1962,  pp. 196-210). 

Applying  the  same  kind  of  approximation  as  in  section 
5.2.3,  equations  (8.1)  yield: 
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where  N  is  the  total  number  of  elementary  areas  in  which 
the  sphere  R  is  subdivided.  The  index  i  refers  to  the 
center  of  the  block  ASj  (this  is  an  approximation  too, 
since  the  point  to  evaluate  the  kernels  in  (8,1)  is  unknown) . 
The  aximuth  ai  ,  the  spherical  distance  ,  and  the  deriva¬ 
tives  of  the  Stokes'  function  with  respect  to  r  ,  and  , 
are  given  by  equations  (6.26)  through  (5.31).  The  differences 
between  equations  (8.2)  of  the  classical  approach,  and  equa¬ 
tions  (5.25)  of  the  mean-value  approach  (following  Bjerhamroar) , 
are: 

(a) .  The  geosphere  Rg  in  (5.25),  now  becomes  the  mean- 

earth  sphere  R  ,  and 

(b) .  The  gravity  anomalies  Ag*  in  (5.25)  do  not  have 

any  physical  meaning  (they  just  satisfy  the  surface 
data),  while  the  anomalies  in  (8.2)  are  the  given 
free-air  anomalies. 

Let  us  now  apply  equations  (8.2)  using  the  data  (surface 
gravity  anomalies)  generated  by  the  terrain  models,  similarily 
to  the  simulations  for  the  improved  techniques.  Using  a 
lO^-inclination  conical  model,  and  gravity  anomalies  on 
a  2'  grid  within  an  area  8®x8®,  the  errors  in  the  components 
of  t  at  10  km  altitude,  are  by  a  factor  of  2  larger  than 
those  from  the  Green's  approach  (table  8.1  vs.  table  4.4). 

The  same  difference  was  found  using  data  on  a  1'  grid  within 
a  0.4®x 0.4®  area  (table  8.3  vs.  4.14). 


Comparing  the  errors  from  the  Dirac  approach  at  5  km, 

10  km,  20  km,  and  100  km  altitudes,  to  those  from  the  clas¬ 
sical  approach,  using  identical  model  and  data,  we  found 
them  to  be  2  to  10  times  smaller  (tables  6.5,  6.6,  6.7, 
and  6.8,  versus  tables  8.2,  8.3,  8.4,  and  8.5  respectively). 
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For  an  easier  comparison  among  the  results  from  the 
various  techniques,  we  list  in  table  8.6  the  exact  component 
of  at  selected  points  at  10  km  altitude,  and  the  compon¬ 
ents  computed  from  the  synthetic  data  on  the  surface  of 
a  10°-inclination  conical  model  (within  a  0.8®x  0.8®  area, 
with  2*  grid),  using  both  the  classical  approach  and  the 
improved  techniques.  It  is  obvious  that  the  results  from 
the  improved  techniques  are  closer  to  the  exJ'.ct  values  than 
the  results  from  the  classical  approach,  with  the  exception 
of  the  results  from  collocation  at  points  uway  from  the 
model's  center  (ifo  =  45®,  Xo  =  250®). 

All  of  our  tests  clearly  indicate  that  the  three  tech¬ 
niques  which  take  the  effect  of  the  topography  into  account, 
offer  an  improved  solution  for  the  computation  of  the  compon 
ents  of  "o  in  space,  as  compared  to  the  classical  approach 
where  the  topography  is  neglected. 


Table  8.1;  Comparison  or'  Gravity  Disturbance  Vector  Components  in  Space 
Computed  Rigorously,  and  From  the  Classical  Approach 
(The  Direct  Integration  Method) 

(values  in  mgal ;  height  in  meters) 
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Table  8.6 

Comparison  of  Gravity  Disturbance  Vector 
Components  in  Space  (10  km  altitude). 

Computed  from  the  Improved  Techniques  and 
from  the  Classical  Approach 
(10°-inclination  conical  model;  0?8x0?8  area; 

2'  grid;  values  in  mgal;  height  in  met"rs) 

Green's  Approach: 

I  EXACT:  I  COMPUTED:  I  DIFTEREHCES:  I 


LAT 

LON 

HEICHTI 

RADIAL 

HORIZ  1 

RADIAL 

HORIZ  1 

RADIAL 

HORIZ  1 

43 

230 

0* 

10000. 1 

-36.89 

0.0  1 

-3 1 . 04 

-2.491 

-5.85 

2.491 

43 

230 

2' 

10000. 1 

-34.66 

-4.881  -28.82 

-5.201 

-6.  14 

0.321 

43 

230 

4* 

10000. 1 

-29.48 

-8.211  -24.25 

-7.281 

-5.23 

-0.931 

43 

230 

6' 

tOOOO. 1 

-23.43 

-9.66 1 

'19.56 

-8.251 

-3.87 

-1.411 

45 

230 

8' 

10000. 1 

-17.93 

-9.781  -15.42 

-8.361 

-2.51 

-1.421 

43 

230 

10' 

10000. 1 

-13.31 

-9.171  -11.94 

-7.881 

-1.57 

-1.281 

43 

230 

12' 

10000. 1 

-10.  13 

-8.241 

-9.33 

-7.201 

-0.82 

-1.041 

Dirac  Approach: 


EXACT: 


I  COMPUTED: 


LAT 

LON 

HEICHTI 

RADIAL 

s 

N 

1 

HORIZ  1 

RADIAL 

43 

230 

0* 

10000.  1 

-36.89 

0.0  1  -35.07 

-0.01 1 

.82 

43 

250 

2' 

10000.  1 

-34.66 

-4.881  -33. 11 

-4.731 

•  1 

.55 

43 

230 

4' 

10000.  1 

-29 . 48 

-8.211  -28.10 

-8.001 

•  1 

.39 

43 

230 

6’ 

10000.  1 

-23.43 

-9.661  -22.18 

-9.491 

.  1 

.25 

43 

250 

8' 

10000. 1 

-17.93 

-9.781  -16.74 

-9.631 

»  1 

. 

45 

230 

10' 

10000. 1 

-13.51 

-9.  171  -  12.30 

-8.981 

•  ) 

43 

230 

12’ 

10000. 1 

-10. 15 

-8.241  -8.96 

-8.01 1 

. 

e.ei  I 
-0. 181 
-•.2J I 
I7l 
-t».  I5» 
-0.  l'»', 
-0.231 


LAT 


Collocation  Approach; 


„  exact:  I  COMPUTED:  1  DlFPEnERCES:  I 

LOW  HEICHTI  RADIAL  HORIZ  I  RADIAL  HORIZ  I  RADIAL  HORIZ  I 


43  230  0’ 

43  230  2’ 

43  230  4’ 

43  230  0 ' 

4.'5  230  8’ 

45  230  10’ 

43  230  12' 


10000. I 
10000. I 
10000. I 
10000. I 
10000. I 
10000.  I 
10000.  I 


-36 

.89 

0. 

.0  1 

-44, 

.99 

-0. 

.021 

8. 

10 

-34, 

.66 

-4, 

.881 

-43, 

.  17 

-4, 

.  101 

8. 

51 

-29. 

.48 

-8. 

.21 1 

-38, 

.47 

-6, 

.901 

8. 

99 

-23. 

.43 

-9, 

.661 

-32. 

69 

-8, 

.001 

9. 

26 

-17. 

.93 

-9. 

781 

-27. 

.45 

-7, 

.761 

9. 

53 

-13. 

SI 

-9. 

171 

-23. 

23 

-6. 

,71  1 

9. 

72 

-10. 

.  15 

-8. 

.241 

-20. 

13 

-5, 

.371 

9. 

99 

0.021 
-0.781 
-1.31  I 
-1.681 
-2.01  I 
-2.461 
-2.871 


LAT 


Classical  Approach; 


I  EXACT:  I  COMPUTED:  I  DIFFERERCES:  I 

LOR  HEICHTI  RADIAL  HORIZ  I  RADIAL  HORIZ  I  RADIAL  HORIZ  I 


43 

250 

O' 

10000. 1 

-36.89 

0.0  1  -22. 17 

-0.041 

-14.73 

0.041 

45 

250 

2' 

10000. 1 

-34.66 

-4.881  -21.39 

-2.421 

-13.28 

-2.461 

45 

230 

4’ 

10000. 1 

-29.48 

-8.211  -19.24 

-4.361 

-10.24 

-3.851 

45 

230 

6’ 

10000. 1 

-23.43 

-9.661  -16.40 

-5.601 

-7.03 

-4.061 

45 

230 

8' 

10000. 1 

-17.93 

-9.781  -13.32 

-6.131 

-4.61 

-3.641 

43 

250 

10’ 

10000. 1 

-13.51 

-9. 171  -10.56 

-6.151 

-2.95 

-3.021 

43 

250 

12' 

10000.  1 

-10.  15 

-8.241  -8.22 

-5.841 

-1.92 

-2.401 

Chapter  9 


SOME  TOPICS  OF  SPECIAL  INTEREST  RELATED  TO  THE 
COMPUTATION  OF  THE  GRAVITY  DISTURBANCE  VECTOR  IN  SPACE 

During  our  investigations  on  the  accuracy  of  the  three 
improved  techniques  for  the  computation  of  the  gravity  vector 
in  space  considering  the  topography  of  the  earth,  some  ques¬ 
tions  related  to  the  applicability  of  these  anoroaches  with 
real  data  were  raised. 

More  specifically,  the  use  of  the  Dirac  approach  with 
data  having  large  spacing  was  questioned,  and  in  section 
9.1  below  some  tests  are  described  that  confirm  this  draw¬ 
back  of  the  method.  In  addition,  the  effect  of  the  truncation, 
i.e.  the  neglect  of  the  information  from  the  remote  zones 
is  discussed  (section  9.2,  and  9.3).  Despite  the  fact  that 
the  effect  of  the  topography  is  not  considered  in  our  compu¬ 
tations  for  the  truncation  effects,  we  think  that  the  results 
in  these  two  sections  will  help  us  in  judging  the  applica¬ 
bility  of  the  Dirac  approu-'.h  with  real  data,  using  a  higher- 
degree  reference  field  as  ^.pposed  to  an  ellipsoidal  one. 


9.1  Some  Tests  with  Real  Data  and  the  Dirac  Approach,  for 

Gravity  Anomaly  Computations  on  the  Surface  of  the  Earth 

In  order  to  test  how  well  the  Dirac  (iterative)  approach 
converges  in  a  real-world  application  with  point  gravity 
anomalies  on  the  surface  of  the  earth,  irregularly  distributed, 
we  decided  to  use  the  data  set  described  in  (Sjoberg,  1978, 
p.64),  and  to  compare  our  results  with  his.  This  data  set 
consists  of  87  point  free-air  gravity  anomalies  in  the  Manitoba 
area,  in  Canada,  with  a  mean  spacing  O'*. 5,  in  an  area  2°.5x6 
and  at  a  mean  elevation  400  meters  approximately.  Note 
that  the  RMS  value  of  these  gravity  anomalies  in  14.00  mgal. 
Table  9.1  contains  the  location,  the  elevation,  and  the 
gravity  anomaly  for  each  one  of  these  points. 

From  these  87  anomalies,  the  gravity  spikes  are  iter¬ 
atively  computed  on  the  geosphere,  whose  depth  from  the 
mean  sphere  R  was  selected  -  for  test  purposes  -  to  be 
0,  10,  20,  and  30  km.  The  Gauss-Seidel  iterative  method 
of  section  5.3.1  was  used  (equation  (5.37)),  using  (5.38) 
as  the  convergence  criterion.  The  RMS  difference  between 
the  given  anomalies,  and  those  computed  from  the  87  spikes 
at  the  same  pointsis  given  in  table  9.3.  The  reason  why 
these  RMS  differences  are  not  exactly  zero,  is  due  to  the 
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fact  that  the  spikes  have  been  computed  iteratively,  and 
not  directly.  In  any  event,  these  differences  are  very 
small,  showing  that  within  three  iterations  the  solution 
for  Ag®  converges. 

Table  9.2  is  of  the  same  nature  as  table  9.1,  but  it 
contains  information  about  an  additional  set  of  50  anomalies 
in  the  same  area,  distributed  within  the  87  anomalies  of 
the  first  set.  The  RMS  value  of  the  50  anraalies,  is  13.5226 
mgal.  The  RMS  differences  between  these  anomalies,  and 
those  computed  from  the  spikes  are  also  shown  in  table  9.3. 

As  the  reader  can  see,  these  RMS  differences  are  of  the 
same  magnitude  as  the  RMS  value  of  the  50  anomalies  itself. 
Sjttberg  (ibid,  p.68)  analyzing  the  same  data  set  with  the 
same  method,  found  similar  results,  but  he  did  not  make 
any  comment  about  the  fact  that  under  the  above  circumstances, 
the  resulting  error  of  the  prediction  is  100%! 

At  first,  if  we  use  (5.47)  to  compute  the  radius  of 
the  smallest  sphere  which  will  guarantee  the  convergence 
of  the  iterative  solution  for  iJ'min=‘0°-5  ,  and  h^ax  “ 
meters,  we  find  that  the  depth  of  the  geosphere  is  about 
28.5  km.  This  is  in  a  very  good  agreement  with  the  optimum 
depth  that  Sjoberg  found  for  the  same  data  set  in  applying 
the  Dirac  approach  (see  his  figure  6,  page  67),  which  is 
approximately  30  km. 

One  reason  for  the  RMS  error  to  be  equal  to  the  RMS 
value  of  the  predicted  quantity,  is  that  the  computed  quantity 
is  orders  of  magnitude  smaller  than  its  true  value  (the 
given  surface  value).  In  order  to  verify  this,  we  computed 
the  gravity  anomalies  at  the  surface  from  the  87  spikes, 
along  a  profile  between  two  of  the  given  data  points.  The 
coordinates,  the  elevation,  and  the  computed  anomalies  at 
selected  points  along  this  profile  are  given  in  table  9.4, 
along  with  the  true  values  at  the  end  points.  We  see  that 
the  magnitude  of  the  gravity  anomaly  is  dramatically  reduced 
as  we  go  away  from  either  one  of  the  end  points.  In  the 
middle  of  the  profile,  the  computed  gravity  anomaly  is  prac¬ 
tically  identical  to  the  given  value. 

Our  interpretation  to  this  -  originally  surprising  - 
finding  is  the  following.  As  the  gravity  anomalies  on  the 
geosphere  are  all  zero  except  at  the  carrier  points,  the 
gravity  anomalies  on  the  surface  of  the  earth  will  also 
retain  the  same  gravity  characteristics  of  the  field  that 
generates  them.  In  other  words,  the  spikes  on  the  geosphere 
fit  only  the  surface  data  from  which  they  have  been  com¬ 
puted,  and  that  cannot  be  used  for  the  computation  of  a 
gravity  anomaly  away  from  these  points,  unless  the  surface 
data  is  very  dense  to  permit  an  accurate  prediction.  This 
also  explains  the  large  errors  that  we  found  from  the  simula¬ 
tion  test  with  5 '-grid  (table  6.3),  as  opposed  to  the  errors 
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Table  9.2 


Fvee-Air  Gravity  Anomalies  in  Manitoba,  Canada 
(Test  Data  Set) 

(From:  Sjoberg  (1978)) 

ST  •  LATITUDE  LOHCITUDE  HEIGHT  ANOIl*.i.Y 
(OECllEES)  (METERS)  ( MCAL) 

“>9999 

l<}009  50.39449  91.14833  306.181  9.76 

1601 1  30.74666  90.96165  379.470  10.63 

10499  30.79678  92.06719  423.501  0.82 

16042  30.95090  91.75833  402.031  -9.20 

16012  50.88303  91.01666  37B.257  1.71 

10052  30.59193  92.61716  420.834  12.66 

10047  50.22284  92.04720  350.445  -3.57 

10030  30.21497  92.37631  357.035  -3.37 

15372  50.49666  91.87166  337.223  12.68 

273  50.40833  91.50833  370.027  3.99 

15370  50.75499  91.88333  391.973  7.53 

16C47  30.68333  91.41998  302.219  13.34 

10*96  50.80524  92.14000  396.545  -2.10 

10081  50.21280  93.23965  364.845  7.53 

623  50.30666  93.17999  361.493  4.41 

16002  31.12833  90.06833  380.390  -21.34 
16050  50.29490  91.30990  373.605  -14.60 
15331  30.74500  90.75665  390.449  7.03 

15334  50.42332  90.71666  391.050  2.09 

16056  50.183.33  90.68032  404.774  -15.16 
15304  50.12999  90.23999  417.576  -3.35 

9164  40.91499  90.60001  457.809  -11.23 
13576  49.02196  91.96101  425.501  -21.95 
13715  49.84846  90.40294  442.265  -2.00 

15723  49.74696  90.75410  435.359  -5.18 

15675  49.06776  92.09227  382.524  -23.51 
10019  49.04630  92.39220  366.674  -10.90 
15070  49.92999  91.38333  308.925  -29.92 
10203  49.76701  94.07744  339.969  2.11 

11433  49.62477  94.02734  373.990  -10.26 
10222  49.59979  94.35619  323.008  5.44 

3710  49.43166  96.27499  357.225  5.27 

5546  49.72333  95.24666  330.328  16.62 

10190  49.62842  95.49942  318.211  -1.36 

5531  49.71666  94.93666  359.664  10.92 

5076  49.71500  94.80666  345.643  14.12 

10070  49.90807  93.14622  372.770  -8.51 

10104  49.60629  93.07337  379.171  -10.07 
10098  49.30528  93.51170  330.023  6.60 

3004  49.81332  92.97501  355.092  -10.10 
10014  49.62167  92.44067  305.572  10.04 

12002  49.49670  92.69481  403.555  11.76 

10007  49.22151  92.46306  405.384  -23.49 
213  49.14166  92.70332  388.620  -22.67 
13723  49.70807  91.09705  435.559  -0.48 

15736  49.66466  91.81667  419.405  -24.77 
13738  49.24759  91.46304  445.617  -13.99 
10249  48.66917  93.27043  337.710  -6.09 

13652  48.56082  90.73390  447.751  -23.06 
15740  48.83151  90.96706  445.922  -10.72 
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from  a  l*-grid  (table  6.4).  These  results  indicate  that 
the  Dirac  approach  should  not  be  used  for  the  interpolation 
of  gravity  anomalies  on  the  surface  of  the  earth,  unless 
the  data  is  sufficiently  dense  to  ensure  an  accurate  predic¬ 
tion.  Note  that  the  errors  in  tables  6.5  through  6.8  are 
much  smaller,  because  of  the  very  find  (1*)  grid  being  used. 


Table  9.3 


Comparison  of  Input  with  Computed  Anomalies,  at 
Various  Depths  of  the  Geosphere  (Dirac  approach). 

Values  in  mgal. 


Depth 

R-Rb 

(km) 

#  of 
Iterat . 

RMS  Diff. 

( input  - 
computed) 

87  points 

RMS  Diff. 
( input  - 
computed) 
50  points 

0 

3 

.52  X  10-® 

13.5266 

10 

3 

.54  X  10~® 

13.5266 

20 

3 

.49  X  10-® 

13.5266 

30 

3 

.61  X  10"« 

13.5266 

Note:  The  iteration  is  terminated  when  two  successive  values 
for  a  spike  have  a  12-digit  agreement. 


Table  9.4 


Gravity  Anomalies  Along  a  Profile 
longitude 


Computation  of  the 


Station 

latitude 

11423 

49‘ 

’01' 

'13V3 

49 

01 

13.4 

49 

01 

12.0 

49 

01 

48.0 

49 

01 

24.0 

49 

06 

00 

49 

12 

00 

49 

18 

00 

49 

24 

00 

49 

30 

00 

49 

31 

48. 

49 

31 

58.8 

49 

31 

59.9 

10013 

49 

32 

00.2 

92‘ 

*58' 

'  56 VO 

408.432 

92 

58 

56.0 

407 

92 

58 

55.2 

405 

92 

58 

48.0 

404 

92 

54 

00 

402.5 

92 

48 

00 

401 

92 

42 

00 

399 

92 

36 

00 

398 

92 

30 

00 

397 

92 

24 

00 

395 

92 

24 

36. 

393 

92. 

24 

32.4 

389 

92 

24 

35.3 

286 

92 

24 

35.6 

383.743 

anomaly 

(mgal)  (mgal) 
21.05  21.0499996 

21.18 
21.01 
1.10 
0.0065 
0.00017 
-0.00011 
-0.000072 
0.000032 
0.0081 
2.45 
5.71 
6.08 

6.16  6.1600003 
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9.2  Truncation  and  Discretization  Errors 


In  the  preceding  chapters  we  presented  three  methods  for 
the  computation  of  the  gravity  vector  in  space,  from  surface 
data,  without  neglecting  the  topography  of  the  earth.  The 
Green's  method  requires  as  data  mean  gravity  disturbances 
on  the  earth's  surface,  while  the  Dirac  and  t’.e  collocation 
approaches  require  point  gravity  anomalies  a’;  discrete  surface 
stations.  From  the  practical  point  of  view,  none  of  these 
methods  can  operate  on  a  global  data  set,  due  to  the  limita¬ 
tions  in  computer  speed  and  core.  One  possible  solution 
to  this  problem  is  to  use  data  in  a  spherical  cap  -  whose 
center  is  at  the  projection  of  the  space  point  on  the  surface 
of  the  earth  -,  and  to  account  fo”  the  information  from 
the  remote  zones  through  a  set  of  potential  coefficients. 
Clearly,  if  only  the  data  in  the  cap  is  used,  the  resulting 
truncation  error  will  be  smaller  (in  general),  as  the  size 
of  the  cap  becomes  larger.  Therefore,  it  is  very  important 
to  know  the  minimum  cap  size  that  yields  a  certain  level 
of  accuracy  in  computing  a  particular  quantity  (the  components 
of  ^  in  our  case). 

In  addition  to  the  truncation  effect,  the  computed 
quantity  is  also  affected  by  the  fact  that  we  are  dealing 
with  discrete  data.  This  type  of  error  is  a  function  of 
the  grid  spacing  (i.e.  the  block  size)  following  the  Green's 
approach,  or  a  function  of  the  distance  between  the  discrete 
data  points  following  the  Dirac  approach,  or  collocation. 

In  general,  the  truncation  error  is  a  function  of  the  extent 
of  the  area,  and  the  discretization  error  is  a  function 
of  the  density  of  the  data. 

In  our  simulations,  the  synthetic  data  on  the  surface 
of  the  model  is  assumed  to  be  errorless  as  being  rigorously 
computed  from  the  disturbing  masses.  Therefore,  the  tables 
in  the  preceding  chapters  give  the  total  error,  which 

is  caused  by  the  combined  effect  of: 

(a) .  Truncation  and  Discretization  (Dirac  approach). 

(b) .  Truncation  and  Discretization,  missing  second-order 

variations  of  the  topography,  and  center-point  evalua¬ 
tions  of  the  kernels  (Green's  approach). 

For  the  simulation  tests  with  the  Green's  approach, 
we  have  used  relatively  small  working  areas  (2'’x2®  to  8°  8®), 
the  reason  being  that  away  from  the  model's  disturbing  masses 
the  value  of  the  gravity  disturbance  diminishes  very  rapidly 
(figure  3.2  through  3.5).  In  other  words,  the  model  creates 
a  very  local  disturbing  field  and  attempts  to  include  more 
and  more  data  by  expanding  the  size  of  the  working  area 
to  8“x  8®  did  not  improve  the  results.  In  addition.,  grid 
intervals  as  small  as  0.5'  in  the  inner  zone  did  not  result 
in  substantial  reduction  of  the  errors  of  .  These  tests 


118. 

indicate  that  the  errors  from  the  Green's  approach  are  mainly 
due  to  the  missing  second-order  variations  of  the  topography, 
and  to  the  effect  of  the  center-point  evaluation  of  the 
kernals.  These  two  approximations  are  the  dominant  sources 
of  errors  in  the  Green's  approach. 

For  the  simulations  with  the  Dirac  approach  we  have 
used  even  smaller  working  areas  (0.4“x0.4®to  2®x2®).  From 
the  results  in  tables  6.3  and  6.4  we  can  see  that  the  errors 
are  now  dramatically  reduced  by  using  five  times  smaller 
grid,  even  if  the  area  extent  becomes  five  times  smaller. 
These  results  indicate  that  in  our  simulations  for  the  Dirac 
approach,  the  discretization  errors  dominate  the  truncation 
errors.  This  is  due  to  the  very  local  characteristics  of 
themodels,  and  it  should  not  be  generalized  for  an  applica¬ 
tion  with  real  data. 


9.3  The  Relationship  between  the  Truncation  Angle,  and 

the  Altitude  of  the  Space  Point 

It  has  been  found  (cf.  Hirvonen  and  Moritz,  1963,  p.68), 
that  in  applications  of  the  classical  approach  for  the  com¬ 
putation  of  the  radial  component  of  t  in  space,  it  is  suffic¬ 
ient  to  extend  the  integration  only  as  far  as  10  times  the 
altitude  of  the  space  point,  in  order  to  ensure  an  error 
smaller  than  lOX.  However,  this  "rule-of-thumb"  is  valid 
only  for  the  Upward  Continuation  Method  (under  a  planar 
approximation),  and  should  not  be  used  for  the  other  two 
techniques  (i.e.  the  Coating  Method,  and  the  Direct  Integra- 
t ion  Method ) . 

From  the  simulations  of  chapter  8  it  was  concluded 
that  the  effect  of  the  topography  is  quite  significant, 
since  the  errors  in  J  are  reduced  by  a  factor  of  2  or 
more,  when  the  topography  is  considered  in  the  improved 
techniques.  The  question  is  now,  how  far  from  the  computation 
point  we  should  extend  the  integration  using  our  improved 
techniques.  In  other  words,  what  is  the  relationship  between 
the  minimum  radius  (((^o)  of  the  cap  inside  of  which  the  data 
is  given,  and  the  altitude  of  the  space  point  (hp),  such 
that  the  truncation  error  in  is  below  a  certain  limit 
(say  ioo;i;)? 

Clearly,  it  is  out  of  the  scope  of  the  present  study 
to  investigate  the  truncation  errors  when  real  data  is  used. 

In  addition,  the  effect  of  the  discretization  is  quite  a 
challenge  to  be  investigated.  A  procedure  for  the  estimation 
of  the  (global  RMS)  truncation  errors  in  computing  the  vector 
t  at  a  space  point  (using  the  Direct  Integration  classical 
approach),  is  described  by  Shepperd  (1979).  We  decided 
to  use  Shepperd 's  procedure  and  computer  programs,  because 
the  topography  away  from  the  computation  point  does  not 
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significantly  contribute  to  the  components  of  t  ,  due  to 
the  diminishing  magnitude  of  the  terms  9S(r,  ,  and 

9S(r,  ii;)/3r  in  equations  (5.48):  Dirac  approach,  and  (8.2): 
classical  approach. 

Shepperd's  equations  for  the  RMS  truncation  error  of 
t  at  altitude  ire: 


radial  error  component:  e(6e_)*r  T  (Ag)  1^  (9.1) 

’’  n»0  ^ 

horizontal  error  component:  e  (6  g.  )»[f^  ("  gr  )‘ +  (9.2) 

n  y  ^ 

•1  ■ 


n=l 


where  On(Ag)  are  the  anomaly  degree  variance.s: 

n 

o^(Ag)  =  C  =  y  (  n-1  V' 


„  =  V  fn-l)*  I  ♦  g*'  ' 

n  nm  nm 

m-0 


(9.3) 


Y  is  the  mean  value  of  the  normal  gravity  (979800  mgal). 
Cj:  ,  §  are  the  fully  normalized  potential  coefficients, 
referred  to  the  same  normal  field  as  the  gravity 
anomalies . 


and  the  truncation  coefficients  for  the  radial  horizontal 
components  of  j  are  defined  as: 


Qn(r,'l/o)  =  R  (^osi*-)  sin^  d^  0.4) 

9n<r,-^o)  =  ^  (9.5) 


where : 


P'(x)  = 


/I  -  x  = 


d  Pn(x) 
dx 


Recursive  relationships  for  the  evaluation  of  these  two 
kinds  of  truncation  coefficients,  starting  from  the  Moloden- 
skii's  coefficients: 


Qjj(r,il^o)  =  S(r,i{>)  Pjj  (cosip)  sinip  d<ti 


(9.6) 


are  given  in  (Shepperd,  1979),  along  with  IDUTHAN  subrouL  ine:i 
for  their  computation.  (These  subroutines  were  converted 
for  double-precision  computations). 


TRUNCfiTION  ERRORS  (MGRL) 

18.  27.  36.  45.  54. 


TRUNCATION  ERRORS  (MGAL) 


TRUNCATION  ANGLE  (DEGREES) 

Radial  (Y),  Horizontal  (X),  and  Total  (+) 
Truncation  Errors  at  Altitude  hp « 10000  meters 
(Maximum  degree:  400,  Ellipsoidal  Reference 
Field) 
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For  the  computation  of  the  anomaly  degree  variances, 
we  can  use  the  Tscherning-Rapp  (1974)  model: 


c  =  A(n-l) 

(n-2)(n+B) 


n  i  3 


A  =  425.28  mgal^ 

B  =  24 


(9.7) 


Using  (9.1)  and  (9.2),  the  total  truncation  error  at 
altitude  will  be: 

total  error  component:  e(<Sg)  =  [e(6g^)^+  e(6gj^)^]^  (9.8) 


In  figures  (9.1),  (9.2),  (9.3),  and  (9.4)  we  give  the 
truncation  errors  for  the  radial,  the  horizontal,  and  the 
total  components  of  ^  at  5  km,  10  km,  100  km,  and  500  km  alt itude 
respectively,  computed  from  equations  (9.1),  (9.2),  and 
(9.8).  The  summations  in  these  equations  were  carried- 
out  from  n  =  2  to  n^ax  ~  ^00  with  C2  =7,5  mgal^  (Tscherning 
and  Rapp,  1974),  i.e.  assuming  an  ellipsoidal  reference 
field.  This  chosen  because  summations  to 

higher  degrees  did  not  yield  significantly  different  results. 
Note  that  Shepperd  (ibid)  performed  his  computations  up 
to  Umax  ®  20  only . 

The  results  plotted  in  these  figures  indicate  that 
in  order  to  maintain  a  truncation  error  smaller  than  lOX 
(with  rspect  to  the  RMS  total  magnitude  of  g  at  the  same 
altitude),  we  should  extend  the  integration  over  a  cap  of 
radius  ip  i  given  in  table  9.5  (second  column)  as  a  funcion 
of  the  altitude  of  the  space  point. 


Table  9.5 

Cap  Radius  for  Truncation  Error  Smaller  than  10%  , 

at  Selected  Altitudes  of  the  Space  Point 

Radius  ipo  : 

hp  Ellipsoidal  Ref.  Field  20th-degree  Ref.  Field 


5 

km 

50° 

3® 

25 

10 

km 

O 

O 

3® 

35 

100 

km 

120® 

5® 

30 

500  km 


155® 


30® 


It  must  be  pointed-out  that  the  fpo  estimates  in  table 
9.5  hold  only  for  the  classical  approach  (chapter  8), 
and  for  the  Dirac  approach  (chapter  5),  both  of  which  are 
based  on  Pizzetti's  generalized  formula  (equations  (5.24), 
and  (8.1)).  As  we  can  see,  by  using  an  ellipsoidal  reference 
field,  even  for  a  relatively  low  altitude  (5  km)  the  inte¬ 
gration  must  be  extended  within  a  50°  cap,  in  order  to  retain 
a  relative  error  below  the  10%  level.  In  order  to  verify 
the  correctness  of  our  computations,  we  computed  (from  Shep- 
perd's  equations)  the  effect  of  truncation  on  Sgj.  ,  and 
iSg.  =  i5g\  =  ,  for  a  point  on  the  surface  of  the  earth. 

We'^used  the  same  anomaly  degree  varip.nces  as  in  (Hirvonen 
and  Moritz,  1963,  p.54): 

C2  =  15;  C3  =43;  =30;  Cs  =  Cs  =  c?  =  Cs  =  25  mgal^ 


and  therefore  the  summations  in  equations  (9.1)  and  (9.2)  were 
taken  to  n^ax  ~  ^  results  are  given  in  table  9.6, 

where  we  also  list  the  truncation  effects  on  6gr  (from 
Hirvonen  and  Moritz,  1963,  table  7.2),  computed  as: 


Q‘  c  ]  = 


(9.9) 


The  truncation  effects  on  6gr  ,  and  5g;«  in  Hirvonen 
and  Moritz  (ibid)  were  incorrectly'^ evaluated  (cf.  Hagiwara, 
1972,  p.457).  These  effects  (on  C  and  n )  should  be  com¬ 
puted  as  (ibid,  p.461): 

1  r  3Ag— 


e(5)  = 

e(n)  = 
where 


(9.10) 


1 

2  ycos^ 


I 

n=2 


^n  ■  "^n  n<n+i)  S(cosiP(,)  (cosi|<o)  sinipo 
and  not  as  in  (Hirvonen  and  Moritz,  ibid,  p.45): 


(9.11) 


e(C) 


:(n )  = 


I  Q. 


n=2 


2  ycosH 


S5 


y  Q 


n=2 


For  reasons  of  symmetry  (ibid,  p.49)  we  can  take 
c(5)  =  e(n)  and  using  "'.10)  we  arrive  to  the  following 
equations  which  give  th^  truncation  effects  on  the  components 


iSg 


<5gi 


of  t 


f 
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£(6g^) 


00 

=  e<Sgj)  .  (i  I  ncn-l)  Q.  =  c„) 
n*2 


i 


(9.12) 


The  results  in  table  9.6  clearly  indicate  that  the 
truncation  effects  on  the  horizontal  components  of  com¬ 
puted  from  Shepperd's  equations  are  identical  to  those 
computed  from  (9.12)  using  Hagiwara's  QjJ  coefficients. 

In  addition,  the  truncation  effects  e(5gj.)  on  the  radial 
component  of  ^  from  Shepperd's  equation  are  in  a  good 
agreement  with  the  Hirvonen-Moritz  (ibid,  table  7.2)  results. 

We  have  also  computed  the  RMS  magnitude  of  the  compon¬ 
ents  of  ^  using  subroutine  COVAX  (Tscherning,  1976).  The 
quantities  that  we  have  actually  computed  form  COVAX,  are 
the  variances  of  — ^  (in  Eotvos),  C  »  and  n  (in 

arcseconds),  i.e.  the  variances  of  the  radial  and  of  the 
horizontal  components  of  ^  ,  at  selected  altitudes  in  space. 
From  these  variances  we  have  then  computed  the  RMS  magnitude 
of  the  components  of  Z  as: 


RMS(6gp)  =  r(var(-(l/i- )3T/9r))^  •  10~’ 
RMS(6gjj)  =  SggSgj'.ggg  (var(S)  +  var(n))^ 
RMS(6gtotai>  =  [RMS(6g^)'+  RMS(6gj^)" 


These  values  must  be  the  same  with  the  truncation  error 
components  (9.1),  (9.2),  and  (9.8),  for  =  0®  ,  provided 
that  the  same  anomaly  degree  variance  model  is  used  (equation 
(9.7)  in  our  case).  The  radial,  the  horizontal,  and  the 
total  components  of  t  from  COVAX,  and  from  Shepperd's 
procedure  (with  nmax^^OO)  are  given  in  table  9.7.  As 
we  can  see,  the  agreement  between  the  results  from  the  two 
procedures  is  very  remarkable. 


Table  9.6 


Comparison  of  the  Truncation  Effects  on  ,  Computed 
from  Shepperd's  Equations,  and  from  Two  Other  Sources 
('^max^S  ,  values  in  mgal,  altitude:  0  meters). 


Cap 

(  +  ) 

( ++ ) 

Shepperd ' 

s  Equations 

radius 

c(<5gj.) 

e(6g^)=  c(^g^) 

e(<^gr> 

c  (f  g-  )  =  c{Sf 

0° 

— 

15.00 

24.11 

15.00 

10“ 

7.5 

10.74 

7.73 

10.74 

O 

O 

5.4 

7.33 

5.89 

7.33 

30“ 

4.4 

5.14 

5.38 

5.14 

o 

O 

4.2 

3.67 

5.59 

3.67 

50“ 

4.2 

2.74 

5.79 

2.74 

60“ 

4.3 

2.29 

5.70 

2.29 

o 

o 

4.1 

2.11 

5.35 

2.11 

o 

O 

00 

3.8 

2.17 

4.84 

2.17 

90“ 

3.3 

2.33 

4.15 

2.33 

100“ 

2.7 

2.44 

3.42 

2.45 

110“ 

2.3 

2.45 

2.93 

2.45 

120“ 

2.2 

2.22 

2.89 

2.22 

130“ 

2.4 

1.77 

3.14 

1.77 

140“ 

2.4 

1.17 

3.28 

1.17 

150“ 

2.2 

0.61 

2.95 

0.61 

160“ 

1.6 

0.20 

2.05 

0.21 

170“ 

0.6 

0.01 

0.74 

0.02 

180“ 

0. 

0. 

0. 

0. 

(+)  from  Hirvonen  and  Moritz,  1963,  table  7.2,  p.55,  using 
the  Molodensky's  Qjj  coefficients. 

(++)  from  equation  (9.12),  using  Hagiwara's  QJ  coefficients 


Table  9.7 
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Comparison  of  the  RMS  Components  of  1  (from  COVAX), 
with  the  Truncation  Errors  for  g-o  =0® 


( t  rom 

Shepperd’s  equ.-* ; 

( value.s 

L  ion*? )  at  Va  T- 1  r  1  ■ 
uro  in  mgal) 

J  •  ’  l; 

hp 

Tscherning ' s 

COVAX 

Shepperd's  equations 

5gr 

6gh 

"^^total 

figr 

‘^^total 

5  km 

38.78 

37.35 

53.84 

37.58 

36.07 

■‘’.09 

10  km 

35 . 76 

34.29 

49.55 

35.53 

33.97 

49.16 

100  km 

23.60 

22.43 

32.56 

23.67 

21.82 

32.20 

500  km 

13.57 

13.91 

19.44 

13.54 

12.01 

18.10 

The  use  of  a  higher  degree  reference  field: 

It  must  be  pointed-out  that  the  large  cap  sizes  (g^o  ) 
given  in  table  9.5  (column  2),  corresponds  to  the  use  of 
;in  ellipsoidal  reference  field.  Let  us  now  assume  ihat 
a  higher-degree  reference  field  is  available,  defined  by 
a  set  of  potential  coefficients  complete  to  degree  and  order 
Href  •  the  GEM-9,  or  the  RAPP-180  solution  (Rapp,  1980). 

If  such  a  reference  field  is  used  for  the  computation  of 
the  components  of  6  at  a  space  point,  the  co'responding 
truncation  effects  can  be  still  computed  from  equations 
(9.1),  (9.2),  and  (9. S'*,  where  the  summations  are  taken 
from  npef  1  Umax  *  example,  if  n^ef  =  20, 

we  found  that  ^  p  ordr^r  to  maintain  the  truncation  errors 
below  the  lOX-ievei,  the  cap  radii  (i^o)  are  much  smaller 
than  the  radii  using  an  ellipsoioal  reference  field  (table 
9.5,  3rfi  column).  Of  course,  these  radii  have  been  computed 
under  the  assumption  the  reference  field  (the  coeffic¬ 

ient.;  IS  error-free.  Nevertheless,  we  can  conclude  that 
(hf  use  of  a  hii;!' -  r-d -gree  reference  field  can  reduce  dram- 
niically  the  truncation  angles,  and  therefore  :iie  computa¬ 
tional  effort. 


SUMMARY,  CONCLUSIONS.  AND  RECOMMEND AT IONS 


Three  approaches  have  been  investigated  xor  the  computa¬ 
tion  of  the  components  of  the  gravity  vector  in  space  from 
surface  data,  without  neglecting  the  topography  of  the  earth: 

1.  A  numerical  integration  approach,  based  on  the  applica¬ 
tion  of  Green's  third  identity  (chapter  2). 

2.  The  Dirac  approach,  following  the  Bjerhammar's  discrete 
solution  to  the  geodetic  B.V.F.  (chapter  5). 

3.  The  Least-Squares  Collocation  approach  (chapter  7). 

In  order  to  avoid  the  errors  which  exist  in  real  data, 
this  work  has  been  a  simulation  study,  using  as  terrain 
a  conical  and  a  spherical  model  (chapter  3).  Seven  para¬ 
meters  are  needed  to  define  the  geometric  and  the  dynamic 
characteristics  of  these  models,  and  the  two  point  masses 
located  beneath  their  surface  on  their  axes,  generate  the 
synthetic  data  for  the  application  of  the  simulations.  The 
exact  gravity  disturbance  vector  components  which  are  computed 
from  the  model ,  are  then  compared  to  the  corresponding  compon¬ 
ents  which  are  evaluated  from  each  one  of  the  three  approaches, 
in  terms  of  their  relative  percentage  differences  (chapters 

4.  6,  and  7).  From  the  simulations  which  were  performed, 
the  following  conclusions  can  be  made. 

(a) .  The  errors  from  the  Green's  approach  become  smaller 
as  the  altitude  of  the  space  point  increases.  They  can 
be  as  small  as  1%  for  a  very  smooth  topography  (the  case 

of  the  spherical  model),  or  for  a  very  large  but  local  topo¬ 
graphic  feature  (the  case  of  the  40°-inclination  cone,  4.1 
km  high).  The  errors  can  be  as  large  as  25%  in  certain 
cases  (right  above  a  20'’-inclination  cone),  but  they  decrease 
at  space  points  away  from  the  model.  These  errors  are  due 
to  the  numerical  integration  procedure  (evaluation  of  the 
kernels  at  the  centers  of  the  blocks),  and  to  the  neglect 
of  the  second  -  order  variations  of  the  topographic  surface. 

(b) .  Both,  the  Dirac  and  the  Collocation  approaches,  require 
discrete  data  on  the  physical  surface  of  the  earth,  but 

they  are  very  time-consuming,  especially  when  a  large  amount 
of  data  is  used.  An  iterative  procedure  for  the  analytical 
continuation  of  the  surface  gravity  anomalies  to  the  geosphere 
is  described,  based  on  the  Gauss-Seidel  numerical  method. 
Acceleration  techniques  which  yield  a  much  faster  rate  of 
convergence,  require  an  additional  series  of  iterations 
for  the  estimation  of  the  eigenvalues  of  the  system  of  equations. 
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which  has  to  be  solved  for  the  gravity  anomalies  on  the 
geosphere.  The  iteration  procedure  might  diverge  for  very 
dense  data  at  high  elevations.  A  method  for  the  estimation 
of  the  depth  of  the  geosphere  is  described  thqt  the 

ircration  is  guaranteed  to  converge.  lb?  i--  t 
gravity  vector  from  the  Dirac  approach  now  increase  as  the 
inclination  of  the  model  increases,  and  as  the  altitude 
of  the  space  point  increases,  but  they  are  still  smaller 
than  the  errors  from  the  Green's  approach  using  identical 
models  and  data. 

(c).  The  collocation  approach  requires  the  computation 
of  the  covariances  between  the  predicted  quanti  t.i  Co ,  ■'iid 
the  data.  In  addition,  it  is  a  very  time-consuming  method, 
since  the  inverse  of  the  matrix  of  the  covariances  is  needed. 

On  an  overall  basis,  the  Dirac  approach  seems  to  be  superior 
as  compared  to  collocation. 

None  of  the  three  approaches  ignore  the  topography 
of  the  earth  as  it  happens  with  the  classical  approach  (chapter 
9).  Comparisons  of  the  improved  techniques  versus  a  classical 
approach  (the  Direct  Integration  Method),  using  identical 
models  and  data  arrangement,  indicated  that  the  errors  in  ? 
from  the  improved  techniques  are  at  least  by  a  factor  of 
2  smaller  than  the  errors  from  the  classical  approach.  This 
clearly  indicates  the  significance  of  the  techniqes  described 
in  this  paper. 

From  some  tests  which  were  made  for  gravity  anomaly 
computations  on  the  surface  of  the  earth  u.sing  real  data 
section  9.1),  it  was  concluded  chat  the  Dirac  approach 
requires  very  dense  data  coverage  for  such  a  kind  of  compu¬ 
tations  (say  1'  as  in  ou:  simulations). 

In  order  i'"<^stigate  the  relationship  between  the 
altitude  of  the  space  point,  and  the  truncation  angle  (the 
radius  of  trie  cap  within  which  the  data  is  given,  section 
9.3),  Shepperd's  (1979)  computer  program  has  been  used. 

It  wa.s  found  that  in  order  to  maintain  a  truncation  error 
in  ^  smaller  than  i0%  twith  respect  to  the  RMS  magnitude 
of  Mie  components  .  r  «  using  an  ellipsoidal  reference 
f  cjd),  we  should  extend  the  computations  within  a  50°- 
ap,  or  larger,  depending  on  the  altitude  of  the  space  point. 
However,  by  using  a  higher-degree  reference  field  (defined 
by  a  set  of  potential  coefficients),  the  truncation  angle 
(.an  be  dramatically  reduced.  The  truncation  effects  discussed 
in  sec’  ion  9.3  are  valid  for  the  Dirac  approach,  .ind  for 
the  •  .assical  approach  (the  Direct  Integration  method), 
bor'i  of  which  are  based  on  the  Pizzetti's  formula. 

From  our  experience  with  the  simulations  pcMformcd 
1  i  described  in  this  study,  we  would  recommend  the  Green's 
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approach  for  the  computation  of  the  gravity  vector  at  high 
altitudes  (above  10  km),  and  the  collocation  approach  for 
points  below  the  10  km-level.  The  Dirac  approach  is  question¬ 
able  due  to  the  fact  that  it  requires  very  dense  data  on 
the  surface  of  the  earth.  We  would  also  recommend  the  use 
of  real  data  (on  the  surface  of  the  earth,  and  in  space), 
for  a  more  realistic  comparison  between  measured  and  computed 
components  of  the  gravity  vector  in  space. 
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